Homogeneous Flux Distribution in High-Flux Solar Furnaces

Comparisons between experimental data and ray-tracing simulation results are presented for the high-flux SF60 solar furnace available at the Plataforma Solar de Almeria, Spain, which has an estimated thermal power of 60 kW. Since an important issue in many applications of solar concentrated radiation is to obtain a radiation distribution that is as homogeneous as possible over the central working area, so-called radiation homogenisers were also used but the degree of success achieved is just satisfactory, as the results show. Finally, further modelling studies using ray-tracing simulations aiming to attain a homogenous distribution of flux by means of double reflexion using two paraboloid surfaces are presented.


Introduction
Solar radiation is an endless source of CO 2 emission-free energy at a global level. Solar heat may not only be used for the production of electricity (especially when integrated with thermal energy storage) but it can also be used in industrial processes [1]. Suitable for installation in zones with high direct normal irradiance, high-flux solar furnaces are a particular type of optical system that use concentrated solar radiation for attaining the high temperatures needed for a wide range of applications in the fields of thermochemistry [2][3][4], materials processing [5,6] and characterisation of materials behaviour [7][8][9][10][11] and properties under extreme conditions [12,13]. The work published in 2017 by Levêque et al. [3] provided a detailed inventory of existing high-flux optical systems and reviewed the designs and their characteristics. It is noteworthy that all existing solar furnaces that make use of real sunlight have been designed and constructed individually, because there are several possible optical configurations [14] which must take into account the geographical location [15] and the maximum power required.
Despite its advantages, there is a difficulty inherent in the current use of point-focussing solar concentration technology for certain applications due to the fact that the flux of solar radiation that reaches the target is theoretically non-homogeneous. As shown in our previous work [14], a circle illuminated with a higher concentration in the middle is theoretically obtained when the paraboloid reflection model [14] is applied to the traditional solar concentrators that use a point-focusing solar concentrating panel assembly (which usually consists of a reflective two-stage concentrator, commonly known as a high-flux solar furnace). The measurements of radiation flux were made with a radiometer Vatell Thermogage 1000-4 (sensitivity = 2 mV per W/cm 2 ) with a circular measurement area (diameter = 5/8 in = 15.9 mm). This is an all-metal circular-foil heat flux transducer, also known as a Gardon-type radiometer. The error described by Ballestrıń et al. [21] originating from the difference in spectrum during calibration and the solar spectrum was not considered. The radiometer was mounted in a moving xyz table, permitting to set the measuring position with an accuracy below 1 mm in all three perpendicular directions.
The experimental measurements were compared with data produced by computer ray-tracing simulations, using the mathematical and computing procedures previously described [14]. Essentially, 10 8 -10 9 parallel rays are generated and then reflected, first by the paraboloid concentrator and then by the 45°-tilted mirror, before being collected by detectors that accumulate the number of The measurements of radiation flux were made with a radiometer Vatell Thermogage 1000-4 (sensitivity = 2 mV per W/cm 2 ) with a circular measurement area (diameter = 5/8 in = 15.9 mm). This is an all-metal circular-foil heat flux transducer, also known as a Gardon-type radiometer. The error described by Ballestrín et al. [21] originating from the difference in spectrum during calibration and the solar spectrum was not considered. The radiometer was mounted in a moving xyz table, permitting to set the measuring position with an accuracy below 1 mm in all three perpendicular directions. The experimental measurements were compared with data produced by computer ray-tracing simulations, using the mathematical and computing procedures previously described [14]. Essentially, 10 8 -10 9 parallel rays are generated and then reflected, first by the paraboloid concentrator and then by the 45 • -tilted mirror, before being collected by detectors that accumulate the number of rays reaching each pixel of the detector screen. This number of rays can then be converted to flux values, for example in kW/m 2 , using standard calibration procedures.
An important goal in many applications of solar concentrated radiation is to obtain a radiation distribution that is as homogeneous as possible over the central working area. Therefore, in the second part of this work, we also used two so-called homogenisers, with 4 and 8 lateral mirrors (see Figure 2), to study the effect produced by these devices on the radiation distribution. No ray-tracing simulations were made for these relatively small devices, because the focal differences between experimental and simulation conditions (a focal diameter of 260 mm versus a focal point) would make it impossible to compare experimental and simulated results.
Energies 2020, 13, x FOR PEER REVIEW 3 of 19 rays reaching each pixel of the detector screen. This number of rays can then be converted to flux values, for example in kW/m 2 , using standard calibration procedures. An important goal in many applications of solar concentrated radiation is to obtain a radiation distribution that is as homogeneous as possible over the central working area. Therefore, in the second part of this work, we also used two so-called homogenisers, with 4 and 8 lateral mirrors (see Figure 2), to study the effect produced by these devices on the radiation distribution. No ray-tracing simulations were made for these relatively small devices, because the focal differences between experimental and simulation conditions (a focal diameter of 260 mm versus a focal point) would make it impossible to compare experimental and simulated results. In the third part of this work, we present the simulations supporting a new type of light concentrating system, based on two paraboloid reflectors, producing an output radiation distribution that is essentially flat over a significant working area, even when non-parallel sun rays are considered. To obtain these non-parallel rays, a possible approach is to generate random points over the Sun's hemisphere illuminating the Earth, using a spherical emission model (see reference [14]) to obtain the polar θ angle ( 01 R = random number between 0.0 and 1.0): and then calculating the deviation angle α from the axis Sun-Earth (between 0° and 0.267°), to obtain the wavelength vector (see Figure 3): In the third part of this work, we present the simulations supporting a new type of light concentrating system, based on two paraboloid reflectors, producing an output radiation distribution that is essentially flat over a significant working area, even when non-parallel sun rays are considered. To obtain these non-parallel rays, a possible approach is to generate random points over the Sun's hemisphere illuminating the Earth, using a spherical emission model (see reference [14]) to obtain the polar θ angle (R 01 = random number between 0.0 and 1.0): and then calculating the deviation angle α from the axis Sun-Earth (between 0 • and 0.267 • ), to obtain the wavelength vector (see Figure 3): tan α = Sun radius × sinθ/distance Sun-to-Earth (2) Energies 2020, 13, x FOR PEER REVIEW 4 of 19 Figure 3. Earth-Sun geometry used to calculate the deviation angle α to parallel rays from randomly generated angle θ, in spherical coordinates.
An essentially equivalent, slightly simpler approach is to generate random points over a circle (the solar disk), to obtain a projected radial distance R: and then calculating the deviation angle α as before: Both methods assume that the solar intensity is the same all over the solar disk. However, astronomical data [22] show that the actual radial profile of intensity resembles a half-circumference, with higher intensity in the centre and reducing to zero in the periphery. To emulate this profile, the wavelength vector length was then multiplied by the corresponding scaling factor. This is an alternative to multiplying the random points over the illuminated half sphere by the cosine of θ (see Figure 3), the so-called Lambertian emitter model that we also tested, arriving to similar results.
Throughout our simulations we aimed to keep the models as simple as possible to emphasize the key aspects, without being distracted by secondary effects. Therefore, atmospheric influences on the incident solar radiation, like those studied by Buie et al. [23], were not considered.
Generating random points in a circle illuminating the concentrator, we obtain a starting vector which combined with the wavelength vector can simulate a random light ray coming from the heliostat [14]. Working out the reflections with the two paraboloids (and filtering lost radiation) and accumulating the rays that intersect detectors carefully positioned, allows measurement of the output radiation that leaves the system and prediction of its overall efficiency in real applications.

Solar Furnace Flux Distribution
Experimental measurements of radiation flux were made for the SF60 solar furnace and compared with modelling data, to assess the distribution of radiation as a function of radial and axial position, in the experimental working zone after the 45°-tilted mirror. As shown in Figure 4a, data was collected after the concentrated radiation is reflected by the 45°-tilted mirror, at z = 475 mm (the focal distance) and z = 730 mm below the horizontal line defined by the concentrator axis. Earth-Sun geometry used to calculate the deviation angle α to parallel rays from randomly generated angle θ, in spherical coordinates.
An essentially equivalent, slightly simpler approach is to generate random points over a circle (the solar disk), to obtain a projected radial distance R: and then calculating the deviation angle α as before: Both methods assume that the solar intensity is the same all over the solar disk. However, astronomical data [22] show that the actual radial profile of intensity resembles a half-circumference, with higher intensity in the centre and reducing to zero in the periphery. To emulate this profile, the wavelength vector length was then multiplied by the corresponding scaling factor. This is an alternative to multiplying the random points over the illuminated half sphere by the cosine of θ (see Figure 3), the so-called Lambertian emitter model that we also tested, arriving to similar results.
Throughout our simulations we aimed to keep the models as simple as possible to emphasize the key aspects, without being distracted by secondary effects. Therefore, atmospheric influences on the incident solar radiation, like those studied by Buie et al. [23], were not considered.
Generating random points in a circle illuminating the concentrator, we obtain a starting vector which combined with the wavelength vector can simulate a random light ray coming from the heliostat [14]. Working out the reflections with the two paraboloids (and filtering lost radiation) and accumulating the rays that intersect detectors carefully positioned, allows measurement of the output radiation that leaves the system and prediction of its overall efficiency in real applications.

Solar Furnace Flux Distribution
Experimental measurements of radiation flux were made for the SF60 solar furnace and compared with modelling data, to assess the distribution of radiation as a function of radial and axial position, in the experimental working zone after the 45 • -tilted mirror. As shown in Figure 4a, data was collected after the concentrated radiation is reflected by the 45 • -tilted mirror, at z = 475 mm (the focal distance) and z = 730 mm below the horizontal line defined by the concentrator axis.  These measurements were made with a radiometer with a circular measurement area (diameter = 5/8 in = 15.9 mm). The measurements were made in a 2D array with 7 × 7 = 49 positions, separated by 15 mm in both x, y directions, covering a sampling region of 105 mm × 105 mm. This sampling region was chosen to fit the output dimensions of the 4-side homogeniser (105 mm × 105 mm). The position at 730 mm was chosen to fit the height of the homogeniser (250 mm), thus allowing a gap of 5 mm between the homogeniser output and the sampling plane. The results obtained with homogenisers are reported in the next Section 4. The energy flux measured (in kW/m 2 ) for plane distances 475 and 730 is reported in Figure 5, using a relative colour scaling, going from blue (cold) to red (hot).    These measurements were made with a radiometer with a circular measurement area (diameter = 5/8 in = 15.9 mm). The measurements were made in a 2D array with 7 × 7 = 49 positions, separated by 15 mm in both x, y directions, covering a sampling region of 105 mm × 105 mm. This sampling region was chosen to fit the output dimensions of the 4-side homogeniser (105 mm × 105 mm). The position at 730 mm was chosen to fit the height of the homogeniser (250 mm), thus allowing a gap of 5 mm between the homogeniser output and the sampling plane. The results obtained with homogenisers are reported in the next Section 4. The energy flux measured (in kW/m 2 ) for plane distances 475 and 730 is reported in Figure 5, using a relative colour scaling, going from blue (cold) to red (hot).
(a) (b) Figure 5. Energy flux (in kW/m 2 ) measured at z = 475 mm (a) and z = 730 mm (b) below the concentrator axis, using a relative colour scale, from blue (cold) to red (hot). The x axis is horizontal. The vertical y direction is aligned with the paraboloid axis, pointing above (+y) to the paraboloid concentrator, and below (−y) to the heliostat.
These measurements were made with a radiometer with a circular measurement area (diameter = 5/8 in = 15.9 mm). The measurements were made in a 2D array with 7 × 7 = 49 positions, separated by 15 mm in both x, y directions, covering a sampling region of 105 mm × 105 mm. This sampling region was chosen to fit the output dimensions of the 4-side homogeniser (105 mm × 105 mm). The position at 730 mm was chosen to fit the height of the homogeniser (250 mm), thus allowing a gap of 5 mm between the homogeniser output and the sampling plane. The results obtained with homogenisers are reported in the next Section 4. The energy flux measured (in kW/m 2 ) for plane distances 475 and 730 is reported in Figure 5, using a relative colour scaling, going from blue (cold) to red (hot).
As expected, the flux and the changes observed at 475 mm (in the focal plane) are much more significant than at 730 mm (i.e., 255 mm away of the focal plane). The announced focal diameter for the SF60 solar furnace is 260 mm (for approx. 90% power), which explains the absence of a well-defined, tiny, focal spot. Moreover, the beam seems to be slightly off-vertical: in the plane 475 the maximum flux occurs below the centre, while in the plane 730 mm the maximum is above the centre. Finally, for z = 730 mm, the column at x = −15 mm shows a slight decrease in flux, due to some minor misalignment. These issues are to be expected in a complex optical system such as SF60, that also depends of the varying sun position.
The sum of the power collected in the sampling region of 105 mm × 105 mm (obtained from data at 7 × 7 = 49 positions) is 17.6 kW, far from the 69 kW announced as the total power achieved by the SF60 solar furnace. This difference is due to the outer illuminated region, which is not collected in this tiny 105 mm × 105 mm region.

Vertical Beam Modelling
As it was previously explained, the measurements of radiation flux were made with a radiometer (Vatell Thermogage 1000 −4 ) with a circular measurement area with diameter of 5/8 inch i.e., 15.9 mm. In order to simulate the conditions of our experimental measurements, Figure 6 shows the illuminated region and the 49 positions where data is collected, assuming either a vertical beam (Figure 6a,c) or a rotated beam (Figure 6b,d). The vertical beam modelling is discussed in Section 3.1. whereas the rotated beam modelling is discussed in the Section 3.2. The 49 data collection positions are represented only by the external contour of the circular measurement area to show as much illuminated area as possible. The plane 475 mm below the concentrator axis corresponds to the experimental focal circle of SF60, with a focal diameter of 26 cm (for~90% power).
In contrast, our model assumes the existence of a focal point, as predicted from geometrical optics. To bridge the gap between both types of results, we shifted the 475-modelling plane upwards until the illuminated circle covers completely the 105 mm × 105 mm region experimentally sampled by the 49 radiometer measurements (see Figure 6a). The resulting shift, z = +90 mm, was applied to both measurement planes, so the plane 475 is now 90 mm above the focal plane and the plane 730 is 730 − 475 − 90 = 165 mm below the focal plane. This is by no means the only option available but it is simple and introduces just one additional degree of freedom. Using a larger z shift, in order to get a 26-cm fully illuminated circle (the announced focal circle for SF60), the generated change is so drastic that experimental and modelling conditions are no longer comparable. Another option would be to randomly smear the rays in a 26-cm circle and then apply a Gaussian or Lorentzian factor, but this introduces more mathematical degrees of freedom without providing further insight on the flux distribution.
The results obtained (number of rays crossing each pixel of the detector) were scaled between the maximum (2079 kW/m 2 ) and minimum (693 kW/m 2 ) experimental values, so modelling and experimental numbers can be directly compared. The resulting values for both 475 and 730 planes are reported in Figure 7, using the same colour procedure of Figure 5.   . Energy flux (in kW/m 2 ) calculated for a vertical beam, for planes z = 475 mm (a) and z = 730 mm (b) below the concentrator axis (shifted 90 mm below), using a relative colour scale, from blue (cold) to red (hot). The x axis is horizontal. The vertical y direction is aligned with the paraboloid axis, pointing above (+y) to the paraboloid concentrator, and below (-y) to the heliostat.   . Energy flux (in kW/m 2 ) calculated for a vertical beam, for planes z = 475 mm (a) and z = 730 mm (b) below the concentrator axis (shifted 90 mm below), using a relative colour scale, from blue (cold) to red (hot). The x axis is horizontal. The vertical y direction is aligned with the paraboloid axis, pointing above (+y) to the paraboloid concentrator, and below (-y) to the heliostat. The values calculated for plane 475−90 mm show a decay from the maximum value that is similar to the behaviour measured with the radiometer, a trend that is made even more clear plotting the energy flux as a function of the +y (vertical) coordinate, as shown in Figure 8, for the vertical beam discussed in this Section 3.1. and the rotated beam discussed in the next Section 3.2.
Energies 2020, 13, x FOR PEER REVIEW 8 of 19 The values calculated for plane 475−90 mm show a decay from the maximum value that is similar to the behaviour measured with the radiometer, a trend that is made even more clear plotting the energy flux as a function of the +y (vertical) coordinate, as shown in Figure 8, for the vertical beam discussed in this Section 3.1. and the rotated beam discussed in the next Section 3.2. The values calculated for plane 730−90 mm are much more flat than the experimental values. This was to be expected because the simulation model is much more sensitive to changes in the distance from the focal plane than the real furnace. As reported above, changing from a focal point to a circle covering completely (circumscribing) a 105-mm square, requires only a vertical shift of +90 mm in the simulation model. In the real furnace, these changes are somehow more averaged, and the focal point is in reality a 26-cm circle (and even so, covering only 90% of the radiation). Increasing even more the z shift used here (+90 mm), the distance to the focal point decreases for the plane 730 (the results look better, less flat) but increases for the plane 475 (the results look worse, more flat).

Rotated Beam Modelling
The experimental results for planes 475 and 730 show (see Figure 5) that the maximum is about 30 mm below the central point for plane 475 and 30 mm above the central point for plane 730 indicating that the beam coming from the 45° mirror is not vertical. This horizontal shift of 60 mm corresponds to a vertical distance of 730 − 475 = 255 mm, so an estimate of the rotating angle can be made: angle = arctan (60/255) = 13.2°.
Although the reason for this angle is probably related to the concentrator's or heliostat's rather complex alignment, it is simpler (and perhaps more enlightening) to model this rotation assuming a shift in the tilted flat mirror (see Figure 1) from 45° to 45° + 13.2°/2 = 51.6° (an α shift in the mirror produces a 2α shift in the output).
The area illuminated by the rotated beam is no longer a circle (see Figure 6b,d) and requires a slightly larger shift in the 475 modelling plane, z = +120 mm, to get all 7 × 7 radiometer positions illuminated and to simulate the trends evidenced by the experimental data. Therefore, the former plane 475 is now 120 mm above the focal plane and the former plane 730 is 730 − 475 − 120 = 135 mm below the focal plane (see Figure 4c). Moreover, the xy centre must be moved horizontally, y = −100 The values calculated for plane 730−90 mm are much more flat than the experimental values. This was to be expected because the simulation model is much more sensitive to changes in the distance from the focal plane than the real furnace. As reported above, changing from a focal point to a circle covering completely (circumscribing) a 105-mm square, requires only a vertical shift of +90 mm in the simulation model. In the real furnace, these changes are somehow more averaged, and the focal point is in reality a 26-cm circle (and even so, covering only 90% of the radiation). Increasing even more the z shift used here (+90 mm), the distance to the focal point decreases for the plane 730 (the results look better, less flat) but increases for the plane 475 (the results look worse, more flat).

Rotated Beam Modelling
The experimental results for planes 475 and 730 show (see Figure 5) that the maximum is about 30 mm below the central point for plane 475 and 30 mm above the central point for plane 730 indicating that the beam coming from the 45 • mirror is not vertical. This horizontal shift of 60 mm corresponds to a vertical distance of 730 − 475 = 255 mm, so an estimate of the rotating angle can be made: angle = arctan (60/255) = 13.2 • .
Although the reason for this angle is probably related to the concentrator's or heliostat's rather complex alignment, it is simpler (and perhaps more enlightening) to model this rotation assuming a shift in the tilted flat mirror (see Figure 1) from 45 • to 45 • + 13.2 • /2 = 51.6 • (an α shift in the mirror produces a 2α shift in the output).
The area illuminated by the rotated beam is no longer a circle (see Figure 6b,d) and requires a slightly larger shift in the 475 modelling plane, z = +120 mm, to get all 7 × 7 radiometer positions illuminated and to simulate the trends evidenced by the experimental data. Therefore, the former plane 475 is now 120 mm above the focal plane and the former plane 730 is 730 − 475 − 120 = 135 mm below Energies 2020, 13, 433 9 of 18 the focal plane (see Figure 4c). Moreover, the xy centre must be moved horizontally, y = −100 mm, the same value for both 475 and 730 planes, to account for the vertical misalignment (the reflected light is further away from the concentrator, see Figure 6b,d).
The values obtained for modelling planes 475 and 730, using the same scaling procedures applied for the vertical beam, are reported in Figure 9. The plane 475 shows a maximum 15 mm below the centre (the radiometer indicates a maximum 30 mm below the centre) while the plane 730 shows a maximum 15 mm above the centre (the radiometer indicates a maximum 30 mm above the centre). Interestingly, it can be seen in Figure 6b,d that, when the maximum value is below the centre of the 105 mm × 105 mm measurement region (plane 475), the centre of the illuminated area is above (the point reached by a laser beam aligned with the concentrator axis, shown as the blue line in Figure 4), and vice-versa (plane 730), showing how difficult it can be to align this solar furnace just using standard experimental techniques. These variations can be understood investigating Figure 4c, where the gradient colours used for both detectors are taken from Figure 6b,d, in the same conditions: when the maximum red is on the right side, the blue laser line is on the left side (plane 475) and vice-versa (plane 730). The decay from the maximum values is similar between calculated and experimental values, for both 475 and 730 planes, as can be seen in Figure 9. The maximum value calculated for plane 730 (940 kW/m 2 ) is almost equal to the equivalent experimental value (950 kW/m 2 ). The correspondence between calculated and experimental data can also be seen in Figure 8b, for the rotated beam.
Energies 2020, 13, x FOR PEER REVIEW 9 of 19 mm, the same value for both 475 and 730 planes, to account for the vertical misalignment (the reflected light is further away from the concentrator, see Figure 6b,d).
The values obtained for modelling planes 475 and 730, using the same scaling procedures applied for the vertical beam, are reported in Figure 9. The plane 475 shows a maximum 15 mm below the centre (the radiometer indicates a maximum 30 mm below the centre) while the plane 730 shows a maximum 15 mm above the centre (the radiometer indicates a maximum 30 mm above the centre). Interestingly, it can be seen in Figure 6b,d that, when the maximum value is below the centre of the 105 mm × 105 mm measurement region (plane 475), the centre of the illuminated area is above (the point reached by a laser beam aligned with the concentrator axis, shown as the blue line in Figure 4), and vice-versa (plane 730), showing how difficult it can be to align this solar furnace just using standard experimental techniques. These variations can be understood investigating Figure 4c, where the gradient colours used for both detectors are taken from Figure 6b,d, in the same conditions: when the maximum red is on the right side, the blue laser line is on the left side (plane 475) and vice-versa (plane 730). The decay from the maximum values is similar between calculated and experimental values, for both 475 and 730 planes, as can be seen in Figure 9. The maximum value calculated for plane 730 (940 kW/m 2 ) is almost equal to the equivalent experimental value (950 kW/m 2 ). The correspondence between calculated and experimental data can also be seen in Figure 8b, for the rotated beam.
(a) (b) Figure 9. Energy flux (in kW/m 2 ) calculated for a 13.2° rotated beam, for planes z = 475 mm (a) and z = 730 mm (b) below the concentrator axis (shifted 120 mm below), using a relative colour scale, from blue (cold) to red (hot). The x axis is horizontal. The vertical y direction is aligned with the paraboloid axis, pointing above (+y) to the paraboloid concentrator, and below (−y) to the heliostat.

Homogeniser Flux Distribution
Experimental results such as those presented in the previous section show that even when the system is completely aligned, the flux distribution tends to decrease rapidly from the centre (in the radial direction), and from the focal point (in the axial direction), significantly decreasing the range of applications for this source of energy. A solution often used to overcome this problem [24][25][26] is to use a multi-mirror homogeniser (Figure 2), to try to obtain a more homogeneous distribution for the flux of energy, as a result of the many reflections that take place inside the homogeniser, that tend to mix the radiation beam.
In our case two homogenisers were tested, a 4-side device and an 8-side device. The more complex 8-mirror homogeniser, with more junctions, did not resist the concentrated power, as shown in Figure 10. . Energy flux (in kW/m 2 ) calculated for a 13.2 • rotated beam, for planes z = 475 mm (a) and z = 730 mm (b) below the concentrator axis (shifted 120 mm below), using a relative colour scale, from blue (cold) to red (hot). The x axis is horizontal. The vertical y direction is aligned with the paraboloid axis, pointing above (+y) to the paraboloid concentrator, and below (−y) to the heliostat.

Homogeniser Flux Distribution
Experimental results such as those presented in the previous section show that even when the system is completely aligned, the flux distribution tends to decrease rapidly from the centre (in the radial direction), and from the focal point (in the axial direction), significantly decreasing the range of applications for this source of energy. A solution often used to overcome this problem [24][25][26] is to use a multi-mirror homogeniser (Figure 2), to try to obtain a more homogeneous distribution for the flux of energy, as a result of the many reflections that take place inside the homogeniser, that tend to mix the radiation beam.
In our case two homogenisers were tested, a 4-side device and an 8-side device. The more complex 8-mirror homogeniser, with more junctions, did not resist the concentrated power, as shown in Figure 10. The flux results obtained for the 4-mirror device, measured with the same radiometer and calibration as before, are shown in Figure 11. The total power, flux average, standard deviation and coefficient of variation (CV = standard deviation/average) are presented in Table 1, for the three sets of experimental data: z = 475 mm and z = 730 mm, with and without homogeniser. Consulting this table and comparing Figure 11 with the equivalent radiation maps shown in Figure 5, three observations can be made. The first is: the values of radiation flux achieved with the homogeniser were much higher than before, even at the focal length, z = 475 mm. The total power went from 17.6 kW at z = 457 mm (focal plane) to 24.7 kW at z = 730 mm with the homogeniser. This is surprising and can only be explained assuming that the homogeniser (whose entrance was on the focal plane) was collecting very oblique rays that would be lost otherwise. We note that the focal diameter of the SF60 solar furnace (260 mm at 90% power), is larger than the entrance length (125 mm) of the homogeniser. Moreover, the slight vertical misalignment of the beam commented upon before certainly contributed to this result. Figure 11. Energy flux (in kW/m 2 ) measured for plane z = 730 mm below the concentrator axis, after a 4-side homogeniser positioned above. The relative colour scale goes from blue (cold) to red (hot). The x axis is horizontal. The vertical y direction is aligned with the paraboloid axis, pointing above (+y) to the paraboloid concentrator, and below (−y) to the heliostat.
The second observation is: the homogeniser produced an output map of radiation that is much more homogeneous than before. This can be concluded by visual inspection of the maps shown in The flux results obtained for the 4-mirror device, measured with the same radiometer and calibration as before, are shown in Figure 11. The total power, flux average, standard deviation and coefficient of variation (CV = standard deviation/average) are presented in Table 1, for the three sets of experimental data: z = 475 mm and z = 730 mm, with and without homogeniser. Consulting this table and comparing Figure 11 with the equivalent radiation maps shown in Figure 5, three observations can be made. The first is: the values of radiation flux achieved with the homogeniser were much higher than before, even at the focal length, z = 475 mm. The total power went from 17.6 kW at z = 457 mm (focal plane) to 24.7 kW at z = 730 mm with the homogeniser. This is surprising and can only be explained assuming that the homogeniser (whose entrance was on the focal plane) was collecting very oblique rays that would be lost otherwise. We note that the focal diameter of the SF60 solar furnace (260 mm at 90% power), is larger than the entrance length (125 mm) of the homogeniser. Moreover, the slight vertical misalignment of the beam commented upon before certainly contributed to this result. The flux results obtained for the 4-mirror device, measured with the same radiometer and calibration as before, are shown in Figure 11. The total power, flux average, standard deviation and coefficient of variation (CV = standard deviation/average) are presented in Table 1, for the three sets of experimental data: z = 475 mm and z = 730 mm, with and without homogeniser. Consulting this table and comparing Figure 11 with the equivalent radiation maps shown in Figure 5, three observations can be made. The first is: the values of radiation flux achieved with the homogeniser were much higher than before, even at the focal length, z = 475 mm. The total power went from 17.6 kW at z = 457 mm (focal plane) to 24.7 kW at z = 730 mm with the homogeniser. This is surprising and can only be explained assuming that the homogeniser (whose entrance was on the focal plane) was collecting very oblique rays that would be lost otherwise. We note that the focal diameter of the SF60 solar furnace (260 mm at 90% power), is larger than the entrance length (125 mm) of the homogeniser. Moreover, the slight vertical misalignment of the beam commented upon before certainly contributed to this result. Figure 11. Energy flux (in kW/m 2 ) measured for plane z = 730 mm below the concentrator axis, after a 4-side homogeniser positioned above. The relative colour scale goes from blue (cold) to red (hot). The x axis is horizontal. The vertical y direction is aligned with the paraboloid axis, pointing above (+y) to the paraboloid concentrator, and below (−y) to the heliostat.
The second observation is: the homogeniser produced an output map of radiation that is much more homogeneous than before. This can be concluded by visual inspection of the maps shown in Figure 11. Energy flux (in kW/m 2 ) measured for plane z = 730 mm below the concentrator axis, after a 4-side homogeniser positioned above. The relative colour scale goes from blue (cold) to red (hot). The x axis is horizontal. The vertical y direction is aligned with the paraboloid axis, pointing above (+y) to the paraboloid concentrator, and below (−y) to the heliostat. The second observation is: the homogeniser produced an output map of radiation that is much more homogeneous than before. This can be concluded by visual inspection of the maps shown in Figures 5 and 11, but also comparing the values of coefficient of variation (CV) in Table 1 for the three maps. While the two flux distributions without homogeniser present CV values of 0.16 and 0.15 (see Table 1), CV decreases substantially, to 0.09, when the homogeniser is used. Thus the homogeniser is fulfilling its original role of flattening the input radiation distribution.
Finally, the third observation: the homogeniser seems to be increasing the beam misalignment, as the maximum radiation spot is now close to the top right corner. As the light paths are now more complex, more degrees of freedom were added to the system, resulting in more misalignment possibilities.
As a general conclusion, it seems fair to say that, although homogenisers are a reasonable workaround to solve an important difficulty in high-flux solar furnaces, they are not a robust, solid solution for the problem. It would be highly desirable for the industry to find a better solution.

Double Paraboloid Flux Distribution
We now analyse a solar concentrator formed first by a large paraboloid (designated as P1) to collect the solar radiation and then a small paraboloid (designated as P2) that reflects back the radiation, producing a very homogeneous, concentrated beam, if the working conditions are properly optimised. We consider first that the input rays are completely parallel and then we analyse the modifications produced in the model when real, slightly non-parallel, sun radiation is considered. Figure 12 shows the general layout for the double paraboloid device, including the radiation detected in various significant locations, for simulations with 10 9 input rays. The key condition is: the two paraboloids must have the same focal point. In this example, the focal point is at z = 0, the small paraboloid P2 is at z = 100 mm and the large paraboloid P1 is at z = 1000 mm. For simplicity, the rim angle was set equal to 45 • . First (stage 1 in Figure 12), the parallel rays (with a central hole as paraboloid P2 blocks the radiation) are reflected by the large, concave paraboloid. The reflected radiation becomes more and more concentrated (stage 2 at z = 500 mm and stage 3 at z = 250 mm) until reaching the small, convex paraboloid, before the focal point (layout A). Then a second reflection occurs and the reflected beam (stage 4 at z = 500 mm and stage 5 at z = 1000 mm) is completely parallel and homogeneous, if a perfect geometry is assumed. The output radiation is collected through a circular hole (with the projected area of paraboloid P2) in the large paraboloid P1. The same output beam can be achieved replacing the convex paraboloid P2 by a concave paraboloid, positioned at the same distance after the focal length, as shown in Figure 12 (layout B). Assuming the rays are parallel, layouts A and B produce exactly the same beam output. Layout A is more compact (and thus produces less diverging errors when the geometry and rays are not mathematically perfect), but layout B uses only concave paraboloids (probably easier to fabricate) and gives access to the true focal point. The parallel, homogeneous output beam is not affected by the position and curvature of the two paraboloids (as long as they have the same focal point), only the beam radius and intensity change (decreasing the size of the second paraboloid produces a beam more intense and narrow). This is a really interesting optical design, probably known from the advent of geometrical optics. However, it requires a high accuracy in the construction of the two paraboloids, something that could hardly be achieved in the past due to constraints in the manufacturing processes.

Parallel Rays
In practise, the geometry is never going to be mathematically perfect in a real life device, so it is important to assess how the behaviour of the device is affected by geometrical inaccuracies. In this work we analysed how a small shift between the position of the two paraboloids affects the output beam produced.
To analyse the beam evolution after the second reflection at paraboloid P2, we prepared four detectors, D1 to D4, as shown in Figure 13. These detectors are equally spaced, at positions z1 = 100 + (1/3)1000 = 433 mm, z2 = 100 + (2/3)1000 = 766 mm, z3 = 100 + (3/3)1000 = 1100 mm and z4 = 100 + (4/3)1000 = 1433 mm. Detectors D1 and D2 are positioned before paraboloid P1 while detectors D3 and D4 are collecting data after the rays went through the hole in P1. We will use this setup from now on, to characterise the output beam in various conditions. By default, we will be using layout A (with a convex paraboloid P2 before the focal point). The parallel, homogeneous output beam is not affected by the position and curvature of the two paraboloids (as long as they have the same focal point), only the beam radius and intensity change (decreasing the size of the second paraboloid produces a beam more intense and narrow). This is a really interesting optical design, probably known from the advent of geometrical optics. However, it requires a high accuracy in the construction of the two paraboloids, something that could hardly be achieved in the past due to constraints in the manufacturing processes.
In practise, the geometry is never going to be mathematically perfect in a real life device, so it is important to assess how the behaviour of the device is affected by geometrical inaccuracies. In this work we analysed how a small shift between the position of the two paraboloids affects the output beam produced.
First we analysed what happens when the distance between the two paraboloids is 1 mm longer (for example setting P2 at z = 99 mm, instead of 100 mm). In Figure 14 we show the radial profiles obtained from data collected at detectors D1, D2 (before P1) and D3, D4 (after P1), for correct (in blue) and incorrect (in red) geometries. In this last case, all rays coming from paraboloid P1 are received and reflected at paraboloid P2, without using its outer border region. The output beam is no longer perfectly homogeneous and converges slightly, being progressively more concentrated from the periphery to the centre. In this case there is no additional loss of energy, as all rays go through the hole at paraboloid P1 and reach the detectors D3 and D4. First we analysed what happens when the distance between the two paraboloids is 1 mm longer (for example setting P2 at z = 99 mm, instead of 100 mm). In Figure 14 we show the radial profiles obtained from data collected at detectors D1, D2 (before P1) and D3, D4 (after P1), for correct (in blue) and incorrect (in red) geometries. In this last case, all rays coming from paraboloid P1 are received and reflected at paraboloid P2, without using its outer border region. The output beam is no longer perfectly homogeneous and converges slightly, being progressively more concentrated from the periphery to the centre. In this case there is no additional loss of energy, as all rays go through the hole at paraboloid P1 and reach the detectors D3 and D4. Figure 14. Computer simulations showing radiation profiles collected at detectors D1, D2 (before P1 hole), D3 and D4 (after P1 hole). In blue for the correct geometry; in red for P2 at z = 99.0 mm, so the distance between P1 and P2 is 1 mm longer. The vertical axis represents the simulated radiation flux (resulting from 10 9 input rays). The horizontal axis represents the horizontal radial distance (with the central hole in the middle).  First we analysed what happens when the distance between the two paraboloids is 1 mm longer (for example setting P2 at z = 99 mm, instead of 100 mm). In Figure 14 we show the radial profiles obtained from data collected at detectors D1, D2 (before P1) and D3, D4 (after P1), for correct (in blue) and incorrect (in red) geometries. In this last case, all rays coming from paraboloid P1 are received and reflected at paraboloid P2, without using its outer border region. The output beam is no longer perfectly homogeneous and converges slightly, being progressively more concentrated from the periphery to the centre. In this case there is no additional loss of energy, as all rays go through the hole at paraboloid P1 and reach the detectors D3 and D4. Figure 14. Computer simulations showing radiation profiles collected at detectors D1, D2 (before P1 hole), D3 and D4 (after P1 hole). In blue for the correct geometry; in red for P2 at z = 99.0 mm, so the distance between P1 and P2 is 1 mm longer. The vertical axis represents the simulated radiation flux (resulting from 10 9 input rays). The horizontal axis represents the horizontal radial distance (with the central hole in the middle).

Figure 14.
Computer simulations showing radiation profiles collected at detectors D1, D2 (before P1 hole), D3 and D4 (after P1 hole). In blue for the correct geometry; in red for P2 at z = 99.0 mm, so the distance between P1 and P2 is 1 mm longer. The vertical axis represents the simulated radiation flux (resulting from 10 9 input rays). The horizontal axis represents the horizontal radial distance (with the central hole in the middle).
Then we analysed what happens when the distance between the two paraboloids is 1 mm too short (for example setting P2 at z = 101 mm, instead of 100 mm). The radial profiles obtained are shown in Figure 15, using the same methodology as before. As the paraboloid P2 is too close to P1, part of the rays (1.7%) coming from paraboloid P1 miss the paraboloid P2. Moreover, as the output beam is slightly divergent, a significant proportion of the rays (12.2%) are unable to go through the hole in paraboloid P1, hitting the inner border instead. Increasing 10 mm the radius of the hole is enough to surpass this problem, at the price of decreasing just slightly (0.2%) the collecting area in the P1 paraboloid. part of the rays (1.7%) coming from paraboloid P1 miss the paraboloid P2. Moreover, as the output beam is slightly divergent, a significant proportion of the rays (12.2%) are unable to go through the hole in paraboloid P1, hitting the inner border instead. Increasing 10 mm the radius of the hole is enough to surpass this problem, at the price of decreasing just slightly (0.2%) the collecting area in the P1 paraboloid. Figure 15. Computer simulations showing radiation profiles collected at detectors D1, D2 (before P1 hole), D3 and D4 (after P1 hole). Blue indicates the correct geometry, red indicates a P2 at z = 101.0 mm, so the distance between P1 and P2 is 1 mm too short. The vertical axis represents the simulated radiation flux (resulting from 10 9 input rays). The horizontal axis represents the horizontal radial distance (with the central hole in the middle).

Non-Parallel Rays
In the previous subsection we assumed that solar rays reach the Earth completely parallel. However, due to the large radius of the Sun's disk, solar rays reach the Earth with an inclination angle between 0° and 0.267° (see Section 2 for the mathematical details). This small inclination is enough to change considerably the behaviour of the double paraboloid device. Using exactly the same geometry and modelling conditions of the previous section, the radiation collected by the four detectors D1 to D4 is shown in Figure 16. Clearly, the rays leaving paraboloid P2 are no longer parallel, they are diverging, approximately 2.49°. Due to the mixed rays, the central hole in the radiation distribution (due to the P2 shadow) is rapidly disappearing and is no longer visible in the working zone, at detectors D3 and D4. Figure 15. Computer simulations showing radiation profiles collected at detectors D1, D2 (before P1 hole), D3 and D4 (after P1 hole). Blue indicates the correct geometry, red indicates a P2 at z = 101.0 mm, so the distance between P1 and P2 is 1 mm too short. The vertical axis represents the simulated radiation flux (resulting from 10 9 input rays). The horizontal axis represents the horizontal radial distance (with the central hole in the middle).

Non-Parallel Rays
In the previous subsection we assumed that solar rays reach the Earth completely parallel. However, due to the large radius of the Sun's disk, solar rays reach the Earth with an inclination angle between 0 • and 0.267 • (see Section 2 for the mathematical details). This small inclination is enough to change considerably the behaviour of the double paraboloid device. Using exactly the same geometry and modelling conditions of the previous section, the radiation collected by the four detectors D1 to D4 is shown in Figure 16. Clearly, the rays leaving paraboloid P2 are no longer parallel, they are diverging, approximately 2.49 • . Due to the mixed rays, the central hole in the radiation distribution (due to the P2 shadow) is rapidly disappearing and is no longer visible in the working zone, at detectors D3 and D4.
Energies 2020, 13, x FOR PEER REVIEW 15 of 19 Figure 16. Computer simulations considering non-parallel rays generated according to equations (1) and (2), showing radiation collected at detectors D1, D2 (before P1 hole), D3 and D4 (after P1 hole). The missing external radiation at detector D3 was blocked by paraboloid P1, when going through the hole. The relative colour scale goes from blue (cold) to red (hot).
Together with the 1.0% of rays lost due to the shadow of paraboloid P2 (discussed in the previous subsection), now 0.20% of the rays fall outside paraboloid P1 (rays generated near the border that point outside), 2.5% of the rays fall outside paraboloid P2 (a single focal point no longer exists) and a staggering 21.1% of the rays do not go through the hole in paraboloid P1. This effect is visible at Figure 16: more external (blue) radiation was not collected at detector D3, because it was blocked by paraboloid P1. This loss is particularly important because this is concentrated radiation, potentially damaging the inner border of paraboloid P1, if no action is taken.
To surpass this problem, the hole radius was increased 30 mm. The concentrated radiation hitting the inner border of paraboloid P1 decreases to just 1.2%, which seems acceptable. With this larger hole, 0.9% of the solar radiation goes directly through the hole, without being concentrated by P1 and P2. In total only 5.3% of the incident radiation is lost. Figure 16. Computer simulations considering non-parallel rays generated according to equations (1) and (2), showing radiation collected at detectors D1, D2 (before P1 hole), D3 and D4 (after P1 hole). The missing external radiation at detector D3 was blocked by paraboloid P1, when going through the hole. The relative colour scale goes from blue (cold) to red (hot).
Together with the 1.0% of rays lost due to the shadow of paraboloid P2 (discussed in the previous subsection), now 0.20% of the rays fall outside paraboloid P1 (rays generated near the border that point outside), 2.5% of the rays fall outside paraboloid P2 (a single focal point no longer exists) and a staggering 21.1% of the rays do not go through the hole in paraboloid P1. This effect is visible at Figure 16: more external (blue) radiation was not collected at detector D3, because it was blocked by paraboloid P1. This loss is particularly important because this is concentrated radiation, potentially damaging the inner border of paraboloid P1, if no action is taken.
To surpass this problem, the hole radius was increased 30 mm. The concentrated radiation hitting the inner border of paraboloid P1 decreases to just 1.2%, which seems acceptable. With this larger hole, 0.9% of the solar radiation goes directly through the hole, without being concentrated by P1 and P2. In total only 5.3% of the incident radiation is lost.
In the previous Section 5.1, using parallel rays, we concluded that the device produces exactly the same results using either the layout A (where paraboloid P2 is convex, positioned before the focal point) or B (where paraboloid P2 is concave, positioned after the focal point). This is no longer the case with non-parallel rays, as the slightly more compact layout A produces slightly smaller losses (the rays diverge less). With layout B, the rays falling outside paraboloid P2 increase to 2.8% (from 2.5%) and the rays that do not go through the hole in paraboloid P1 increase to 2.0% (from 1.2%) so the total radiation lost increases to 6.9% (from 5.3%). Unless mentioned otherwise, we will always be using layout A in our analysis, for simplicity.
A more detailed analysis of the radiation reaching the various detectors can be obtained plotting the detected radiation as a function of the radial distance, as shown in Figure 17. Detectors D1 and D2 still show the effect of the P1 shadow, but after going through the hole, at detector D3, the radiation is already perfectly homogeneous over a large central region (with a diameter of 65 mm), with almost the same intensity as obtained in the previous section (with parallel rays). Outside the central region, the intensity decreases almost linearly. The diameter of the central region decreases with distance: it is 65 mm at detector D3 and only 35 mm at detector D4.
Energies 2020, 13, x FOR PEER REVIEW 16 of 19 Figure 17. Computer simulations showing radiation profiles collected at detectors D1, D2 (before P1 hole), D3 and D4 (after P1 hole). Blue indicates parallel rays, red non-parallel, equally distributed rays. Non-parallel rays were generated according to Equations (3) and (4). The vertical axis represents the simulated radiation flux (resulting from 10 9 input rays). The horizontal axis represents the horizontal radial distance (with the central hole in the middle).
The intensity of the concentrated radiation in the central zone (45,000 rays) is now circa 4.4% lower than with parallel rays (due to the larger P1 hole), but all the benefits of the concentrated, homogeneous radiation are still there. In particular, at detector D3, the diameter of the central zone (65 mm) is a very useful distance for many applications. This value is far from the 170 mm measured for parallel rays, but in compensation is a full disk, without the 18 mm central hole. Of course, all the crown region around, with a 60 mm width, is still valuable, concentrated radiation, it is just no longer homogeneous, decreasing linearly from 45,000 rays in the central homogeneous region to just 15,000 rays in the borders.
Clearly, the radiation reaching detector D3 is well suited for many scientific and industry applications, as it is homogeneous over a large central zone and is essentially as concentrated as in Figure 17. Computer simulations showing radiation profiles collected at detectors D1, D2 (before P1 hole), D3 and D4 (after P1 hole). Blue indicates parallel rays, red non-parallel, equally distributed rays. Non-parallel rays were generated according to Equations (3) and (4). The vertical axis represents the simulated radiation flux (resulting from 10 9 input rays). The horizontal axis represents the horizontal radial distance (with the central hole in the middle).
The intensity of the concentrated radiation in the central zone (45,000 rays) is now circa 4.4% lower than with parallel rays (due to the larger P1 hole), but all the benefits of the concentrated, homogeneous radiation are still there. In particular, at detector D3, the diameter of the central zone (65 mm) is a very useful distance for many applications. This value is far from the 170 mm measured for parallel rays, but in compensation is a full disk, without the 18 mm central hole. Of course, all the crown region around, with a 60 mm width, is still valuable, concentrated radiation, it is just no longer homogeneous, decreasing linearly from 45,000 rays in the central homogeneous region to just 15,000 rays in the borders.
Clearly, the radiation reaching detector D3 is well suited for many scientific and industry applications, as it is homogeneous over a large central zone and is essentially as concentrated as in the previous experiment, with parallel rays. Moving further away the working zone, as detector D4 shows, the radiation becomes more concentrated in the homogeneous, central zone, but this zone becomes noticeably smaller.

Calibrated Rays
In the previous subsection we assumed that the incoming radiation intensity was equal all over the Sun disk. Accurate astronomical results actually show that the intensity falls sharply near the border, according to a radial profile from r = −radius to r = +radius that resembles a half-circumference. Applying the resulting correction factor, the device changes slightly in its performance, as shown in Figure 18, for the same modelling conditions. Energies 2020, 13, x FOR PEER REVIEW 17 of 19 Figure 18. Computer simulations showing radiation profiles collected at detectors D1, D2 (before P1 hole), D3 and D4 (after P1 hole). Blue indicates parallel rays, red indicates non-parallel, calibrated rays. Calibrated rays were generated multiplying a scaling factor, as explained in Section 5.3. The vertical axis represents the simulated radiation flux (resulting from 10 9 input rays). The horizontal axis represents the horizontal radial distance (with the central hole in the middle). Now 0.17% of the rays fall outside paraboloid P1 (instead of 0.20%), 2.2% of the rays fall outside paraboloid P2 (instead of 2.5%) and only 0.6% of the rays do not go through the hole in paraboloid P1 (instead of 1.2%). The radiation intensity in the central zone for detector D3 is now circa 6.8% lower than with parallel rays but in compensation the diameter of the homogeneous, concentrated, region increased to 75 mm (instead of 65 mm). At detector D4 the central zone increased to 40 mm (from 35 mm). These changes are relatively small, however, so real data (measured by radiation sensors) might not be too far away from the simulation results obtained here.

Final Remarks
During this work we did not attempt to optimise the device geometry. For example, the influence of the distance between the two paraboloids on the output flux distribution, for the calibrated rays considered in Section 5.3, is not yet studied. Geometrical parameters such as the rim angle, the ratio between the two paraboloid focal distances and the hole aperture, are likely to play a key role in the overall performance. More modelling analysis are recommended, to shed light on these aspects, before attempting to construct this device. Figure 18. Computer simulations showing radiation profiles collected at detectors D1, D2 (before P1 hole), D3 and D4 (after P1 hole). Blue indicates parallel rays, red indicates non-parallel, calibrated rays. Calibrated rays were generated multiplying a scaling factor, as explained in Section 5.3. The vertical axis represents the simulated radiation flux (resulting from 10 9 input rays). The horizontal axis represents the horizontal radial distance (with the central hole in the middle). Now 0.17% of the rays fall outside paraboloid P1 (instead of 0.20%), 2.2% of the rays fall outside paraboloid P2 (instead of 2.5%) and only 0.6% of the rays do not go through the hole in paraboloid P1 (instead of 1.2%). The radiation intensity in the central zone for detector D3 is now circa 6.8% lower than with parallel rays but in compensation the diameter of the homogeneous, concentrated, region increased to 75 mm (instead of 65 mm). At detector D4 the central zone increased to 40 mm (from 35 mm). These changes are relatively small, however, so real data (measured by radiation sensors) might not be too far away from the simulation results obtained here.

Final Remarks
During this work we did not attempt to optimise the device geometry. For example, the influence of the distance between the two paraboloids on the output flux distribution, for the calibrated rays considered in Section 5.3, is not yet studied. Geometrical parameters such as the rim angle, the ratio between the two paraboloid focal distances and the hole aperture, are likely to play a key role in the overall performance. More modelling analysis are recommended, to shed light on these aspects, before attempting to construct this device.

Conclusions
Obtaining homogeneous, concentrated solar radiation, has long been a goal that has not yet been achieved. The comparison between experimental and simulated results for high-flux SF60 solar furnace demonstrated that there is still a long way to go to improve the accuracy and alignment of this large-scale equipment. Even if the output radiation was a perfect Gaussian or Lorentzian function, it would be far from the desired homogeneous output radiation.
Radiation homogenisers have been used as a useful work-around, to mix more the radiation, but the degree of success achieved is just satisfactory, as our results show. A possible solution is to use a device with two paraboloids, with the same focal point, as discussed in this work. With parallel rays, this device works like a fairy tale, but it seems very sensitive to solar ray inclination and the device's manufacturing quality. Manufacturing this device with the required geometrical accuracy requires very precise mechanical engineering, that is likely to increase the costs significantly. More modelling studies are encouraged, to optimise further the geometry, before building a physical and expensive device.