Problem of Heliostat Field Efficiency Calculation Based on Innovative Monte Carlo Simulation

1. Introduction

This study focuses on the efficiency of a heliostat field system for focusing sunlight to a collector, where the average annual thermal power output and optical efficiency are the key metrics for evaluating its performance. These metrics include five core components: shadow shading efficiency, cosine efficiency, atmospheric transmittance, collector truncation efficiency, and mirror reflectance. In this work, a heliostat field located at a specific east longitude, north latitude and altitude with a specific radius is selected as the research object, and the corresponding mathematical model of optical efficiency is established to facilitate the calculation of these five indexes. The calculation of the shadowing and occlusion efficiency is particularly complicated, and we convert the original three-dimensional problem into a two-dimensional planar problem[1] by establishing the transformation matrices between the mirror coordinate system and the mirror field coordinate system, which are processed with the help of the improved Monte Carlo simulation method. Using these transformation matrices for coordinate transformation, ray tracing simulations are performed on the 2D plane, and the shadowing and occlusion efficiencies are successfully calculated. Meanwhile, the parallel computation is accelerated by further improvement of the Monte Carlo algorithm, which improves the computational efficiency of other average optical efficiencies of the fixed-heaven mirror field. In the end, we simplified the calculation of the annual average thermal power output and other related performance indexes by using the five time points of 9:00, 10:30, 12:00, 13:30 and 15:00 on the 21st day of each month as the "annual average indexes".

2. Modeling of Heliostat Field and Efficiency Calculations

2.1. Modeling of Heliostat Field

The distribution of the heliostat field is shown in the figure below, with the center of the absorber tower base as the origin, the due east direction as the positive direction of the $X_g$ axis, the due north direction as the positive direction of the $Y_g$ axis, and perpendicular to the center of the absorber tower base upward as the positive direction of the $Z_g$ axis.

Figure 1. Schematic diagram of the heliostat field coordinate system

Figure 1. Schematic diagram of the heliostat field coordinate system

2.2. Cosine Efficiency

For a single heliostat, the cosine efficiency can reach its maximum value only when the incident ray vector is parallel to the reflected ray vector and the normal vector.The derivation of the cosine efficiency calculation is as follows:

The solar altitude angle ${\alpha }_s$ and solar azimuth angle ${\gamma }_s$ are:

$$\left\{\begin{array}{c}sin{\alpha }_s=cos\delta cos\phi cos\omega +sin\delta sin\phi \hfill \\ cos{\gamma }_s=\frac{sin\delta -sin{\alpha }_{s}sin\phi }{cos{\alpha }_{s}cos\phi }\hfill \\ \omega =\frac{\pi }{12}\left(ST-12\right)\hfill \\ sin\delta =sin\frac{2\pi D}{365}sin\left(\frac{2\pi }{360}23.45\right)\hfill \end{array}\right.$$

(1)

where $\phi$ denotes local latitude, with north latitude being positive; $\omega$ denotes solar time angle, $ST$ denotes local time, $\delta$ denotes solar declination, and D denotes the number of days counted from the vernal equinox as the 0th day.

The position of the sun can be represented by the solar altitude angle ${\alpha }_s$ and the solar azimuth angle ${\gamma }_s$, and for the light incident on the heliostat the vector is represented as $\stackrel{⟶}{G_{r,M}}\left(x_r,y_r,z_r\right)$, with coordinates:

$$\left\{\begin{array}{c}{x}_r=cos{\alpha }_{s}sin{\gamma }_s\hfill \\ y_r=cos{\alpha }_{s}cos{\gamma }_s\hfill \\ z_r=sin{\alpha }_s\hfill \end{array}\right.$$

(2)

As shown in the Figure 2, in the heliostat field , the collector center coordinates for $J\left(0,0,H_T\right)$ , where $H_T$ for the height of the absorber tower, for any heliostat field center $M_i$, its coordinates for $M_i\left(x_{M,i},y_{M,i},z_{M,i}\right)$, where M for the heliostat field center and i denotes the $ith$ heliostat field mirror. Vector $\stackrel{⟶}{G_{f,J}}$ represents the unit vector J pointing to the center of the collector from the center of the heliostat field $M_i$, denoted as:

$$\stackrel{⟶}{G_{f,J}}=\frac{J-M}{\left|J-M\right|}=\frac{\left(-x_{M,i},-y_{M,i},H_T-z_{M,i}\right)}{\sqrt{{x^2}_{M,i}+{y^2}_{M,i}+{\left(H_T-z_{M,i}\right)}^2}}$$

(3)

According to the law of reflection of light, the mirror normal vector $\stackrel{⟶}{G_{n,M}}$ bisects the angle between the incident unit vector $\stackrel{⟶}{G_{r,M}}$ and the reflected unit vector $\stackrel{⟶}{G_{f,J}}$:

$$\stackrel{⟶}{G_{n,M}}=\frac{\stackrel{⟶}{G_{r,M}}-\stackrel{⟶}{G_{f,J}}}{\left|\stackrel{⟶}{G_{r,M}}-\stackrel{⟶}{G_{f,J}}\right|}=\left(o,p,q\right)$$

(4)

Then the heliostat pitch angle A and heliostat azimuth angle B can be derived from the mirror normal vector $\stackrel{⟶}{G_{n,M}}$

$$\left\{\begin{array}{c}cos\left(A\right)=\frac{q}{\sqrt{o^2+p^2+q^2}}\hfill \\ tan\left(B\right)=\frac{p}{o}\hfill \end{array}\right.$$

(5)

As shown in the Figure 2, there is an angle $\theta$ between the incident or reflected ray and the normal, where the cosine value $cos\theta$ is the value of the cosine efficiency of that fixed-sun mirror:

$${\eta }_{cos}=cos\theta =\frac{-\stackrel{⟶}{G_{r,M}}·\stackrel{⟶}{G_{n,M}}}{\left|\stackrel{⟶}{G_{r,M}}\right|\left|\stackrel{⟶}{G_{n,M}}\right|}=-\stackrel{⟶}{G_{r,M}}·\stackrel{⟶}{G_{n,M}}$$

(6)

where $-G_{r,M}$ denotes the opposite direction of incident light:

$$-\stackrel{⟶}{G_{r,M}}=\left(-cos{\alpha }_{s}sin{\gamma }_s,-cos{\alpha }_{s}cos{\gamma }_s,-sin{\alpha }_s\right)$$

(7)

2.3. shadow Shading Efficiency

In heliostat fields, not all incident and reflected sunlight can be perfectly focused on the central receiver, leading to two primary types of optical losses. First, shadow loss occurs when certain heliostats, referred to as obstructing heliostats, block sunlight that should ideally reach a targeted heliostat. Secondly, blocking loss arises when light, reflected off the target heliostat, is intercepted by heliostats positioned in front of it. These phenomena, distinguishing between shadow and blocking losses, are critical to understanding and optimizing the efficiency of heliostat fields.The blocking loss and shading loss of the mirror are shown below:

Figure 3. Schematic diagram of blocking loss and shadow loss for mirrors

Figure 3. Schematic diagram of blocking loss and shadow loss for mirrors

Not all heliostats surrounding a target heliostat are considered obstructing heliostats. According to reference [1], it is sufficient to consider only the heliostats located in the immediate vicinity of the target heliostat. Given that the distance between the centers of adjacent heliostat bases exceeds the width of the mirror surface by 5 meters, the area of concern for obstructing heliostats is defined by a circle centered on the base of the target heliostat with a radius of 11 meters.

As illustrated in Figure 2, we establish a mirror coordinate system for the heliostat, designating the center of the mirror surface, M, as the origin. The $X_m$ axis is defined by a line passing through the origin and parallel to the mirror’s width, while the $Y_m$ axis is defined by a line passing through the origin and perpendicular to the mirror’s width. Additionally, the $Z_m$ axis is defined by a line passing through the origin and perpendicular to the $X_m{Y}_m$ plane, with orientations as depicted in the figure. Our objective is to determine whether a ray of light, encompassing both incident and reflected rays, passing through a point $C_a\left(x_a,y_a,0\right)$ on mirror a, will intersect mirror b. If so, the intersection point on mirror b, denoted as $C_b\left(x_b,y_b,0\right)$, is expressed in the mirror’s coordinate system[1].

For a mirror in a heliosta field, the horizontal axis controls the pitch angle of the mirror and the vertical axis controls the azimuth angle of the mirror, and we can derive the transformation matrix between the mirror coordinate system and the mirror field coordinate system from the pitch angle and the azimuth angle as follows:

$$T=\left(\begin{array}{ccc}-sinA& -sinBcosA& cosBcosA\\ cosA& -sinBsinA& cosBsinA\\ 0& cosB& sinB\end{array}\right)$$

(8)

where A represents the heliostat pitch Angle and B represents the heliostat azimuth Angle, which are given by equation 5.

For a point $C_a\left(x_a,y_a,0\right)$ in mirror A, it is converted from mirror coordinate $C_a\left(x_a,y_a,0\right)$ to coordinate ${C^{\prime }}_a\left({x^{\prime }}_a,{y^{\prime }}_a,0\right)$ in heliostat field coordinate system, the calculation process is as follows:

$${C^{\prime }}_a=\left(\begin{array}{ccc}-sinA& -sinBcosA& cosBcosA\\ cosA& -sinBsinA& cosBsinA\\ 0& cosB& sinB\end{array}\right)·C_a+M_a=\left(\begin{array}{c}{x^{\prime }}_a\\ {y^{\prime }}_a\\ {z^{\prime }}_a\end{array}\right)$$

(9)

where $M_a$ represents the center of mirror a, that is, the coordinates of the origin of the mirror a coordinate system within the field coordinate system, denoted by $\left(x_{M,a},y_{M,a},z_{M,a}\right)$.

The transformation of point ${C^{\prime }}_a\left({x^{\prime }}_a,{y^{\prime }}_a,0\right)$ from the field coordinate system to the coordinate system of mirror b, resulting in ${C^{″}}_a\left({x^{″}}_a,{y^{″}}_a,0\right)$ :

$${C^{″}}_a={\left(\begin{array}{ccc}-sinA& -sinBcosA& cosBcosA\\ cosA& -sinBsinA& cosBsinA\\ 0& cosB& sinB\end{array}\right)}^T·\left({C^{\prime }}_a-M_b\right)=\left(\begin{array}{c}{x^{″}}_a\\ {y^{″}}_a\\ {z^{″}}_a\end{array}\right)$$

(10)

Where $M_b$ represents the center of mirror b, which is the coordinates of the origin of mirror b’s coordinate system within the field coordinate system, denoted by $\left(x_{M,b},y_{M,b},z_{M,b}\right)$.

When a beam of light enters the heliostat coordinate system, its vector representation is denoted as $\stackrel{⟶}{G_M}$ In the field coordinate system, the vector representation is $\stackrel{⟶}{G_P}$. The transformation relationship between them is as follows:

$$\stackrel{⟶}{G_M}={\left(\begin{array}{ccc}-sinA& -sinBcosA& cosBcosA\\ cosA& -sinBsinA& cosBsinA\\ 0& cosB& sinB\end{array}\right)}^T·\stackrel{⟶}{G_P}=\left(e,f,d\right)$$

(11)

Summarizing the above, with the point ${C^{″}}_a\left({x^{″}}_a,{y^{″}}_a,0\right)$ known within the coordinate system of mirror b, along with the light vector $\stackrel{⟶}{G_M}\left(e,f,d\right)$, we can determine the intersection point $C_b\left(x_b,y_b,0\right)$ of the light vector with the plane of mirror b in its coordinate system. The relationship between them is as follows:

$$\left\{\begin{array}{c}\frac{x_b-{x^{″}}_a}{e}=\frac{y_b-{y^{″}}_a}{f}=\frac{0-{z^{″}}_a}{d}\hfill \\ x_b=\frac{{d{x}^{″}}_a-{e{z}^{″}}_a}{d}\hfill \\ y_b=\frac{{d{y}^{″}}_a-{f{z}^{″}}_a}{d}\hfill \end{array}\right.$$

(12)

If $-3\le x_a\le 3$ and $-3\le y_a\le 3$, then point $C_2$ falls within the range of mirror b, indicating it is within the shadow and blocking area.

Therefore, we employ the Monte Carlo simulation algorithm to uniformly generate $n\times n$ points on the rectangular surface of mirror a. Through the calculation process described above, we determine which of these $n\times n$ points fall onto mirror b. Let the number of points falling onto mirror b be denoted by $N_b$. Consequently, the shadow and blocking efficiency, ${\eta }_{sb}$, is defined as follows:

$${\eta }_{sb}=1-\frac{N_b}{n\times n}$$

(13)

2.4. Collector Truncation Efficiency

The collector truncation efficiency refers to the ratio of the solar radiant energy actually captured by the solar collector to the theoretical maximum that could be captured. In practice, there are certain losses in the energy radiated to the collector, such as those due to the precision limitations of the heliostats or the displacement of the focal spot caused by swaying. This concept is essential for understanding the real-world performance of solar thermal systems and identifying areas for technological improvement.

To calculate the collector truncation efficiency, we adopted the method proposed by Collado and Guallar, based on the HFLCAL model [3]:

$${\eta }_{trunc}={\int }_{-R}^R\frac{1}{{\sigma }_{tot}\sqrt{2\pi }}{e}^{\left(-\frac{x^2}{2{{\sigma }_{tot}}^2}\right)}dx{\int }_{\left(-\frac{H_T}{2}-h_0\right)·cos\theta }^{\left(\frac{H_T}{2}-h_0\right)·cos\theta }\frac{1}{{\sigma }_{tot}\sqrt{2\pi }}{e}^{\left(-\frac{y^2}{2{{\sigma }_{tot}}^2}\right)}dy$$

(14)

where R denotes the radius of the collector, and $H_T$ represents the height of the absorption tower, which is the distance from the center of the collector to the ground. $h_0$ refers to the installation height of the heliostat. $\theta$ signifies the angle formed by the incident (or reflected) light rays with the normal. The total impact factor${\sigma }_{tot}$is determined by a combination of factors: ${\sigma }_{sun}$, which represents the solar shape error; ${\sigma }_{bq}$, which accounts for the beam quality error; ${\sigma }_{ast}$, denoting the astigmatism error; and ${\sigma }_{track}$, which stands for the tracking error. The calculation process for these components is as follows:

$${\sigma }_{\mathrm{tot}}=\sqrt{{\sigma }_{\mathrm{sun}}^2+{\sigma }_{\mathrm{bq}}^2+{\sigma }_{\mathrm{ast}}^2+{\sigma }_{\mathrm{track}}^2}$$

(15)

where the values for $\sigma$ are referenced from relevant literature [4]. Specifically, ${\sigma }_{sun}$ is set at 2.51 mrad, ${\sigma }_{track}$ at 0.63 mrad, and ${\sigma }_{bq}={\sigma }_s^2$, where ${\sigma }_s$ denotes the slope error and is valued at 0.94 mrad.

The astigmatism error ${\sigma }_{ast}$ is calculated using the following equation:

$${\sigma }_{\mathrm{ast}}=\frac{\sqrt{0.5\left(H_{\mathrm{t}}^2+W_{\mathrm{s}}^2\right)}}{4D_{\mathrm{rec}}}$$

(16)

where $H_T$ and $W_s$ respectively represent the height and width of the heliostat, while $D_{rec}$ denotes the distance from the center of the heliostat to the center of the collector, which is significantly greater than the dimensions of the heliostat. This substantial difference causes the denominator of the equation, $D_{rec}$, to be much larger than the value of the numerator, leading to the approximation that the value of ${\sigma }_{ast}$ is effectively 0.

2.5. Atmospheric Transmittance and Mirror Reflectance

i.atmospheric transmittance

The atmospheric transmittance ${\eta }_{at}$ in a heliostat field refers to the extent to which sunlight is affected by absorption, scattering, and refraction as it passes through the atmosphere to reach the heliostats. It quantifies the proportion of sunlight that is successfully received by the heliostats [6].

$${\eta }_{\mathrm{at}}=0.99321-0.0001176D_{\mathrm{rec}}+1.97\times {10}^{-8}\times {D^2}_{\mathrm{rec}}\left(D_{\mathrm{rec}}\le 1000\right)$$

(17)

where $D_{rec}$ denotes the distance from the center of the heliostat to the center of the collector.

ii.mirror reflectance

The mirror reflectance, ${\eta }_{ref}$, in a heliostat field denotes the proportion of sunlight that is reflected by the surface of the heliostats upon incidence. It reflects the ability of the heliostats to redirect incoming solar radiation towards the target receiver. This parameter can be considered a constant, with a typical value of ${\eta }_{ref}=0.92$.

2.6. Optical Efficiency and Thermal Power Output

i.optical efficiency

Summarizing the 5 methods of calculating optical efficiency, we can use the following formula to calculate the total optical efficiency [5]:

$$\eta ={\eta }_{cos}{\eta }_{sb}{\eta }_{trunc}{\eta }_{at}{\eta }_{ref}$$

(18)

ii.thermal power output

The output thermal power $E_{fiel}$ of a heliostat field depends on several factors, and it is given by the following equation:

$$E_{\mathrm{field}}=\mathrm{DNI}·\sum _i^N{A}_i{\eta }_i$$

(19)

where $N_{mirror}$ denotes the total number of heliostat mirrors, $DNI$ denotes the normal direct radiation irradiance, ${\eta }_i$ denotes the optical efficiency of the $ith$ mirror, and $A_i$ denotes the light-gathering area of the $ith$ mirror.

Normal direct radiation irradiance is the irradiance of sunlight directed onto the ground in a direction perpendicular to the ground. It represents the amount of solar energy received per unit area, as follows [7]:

$$\left\{\begin{array}{c}\mathrm{DNI}=G_0\left[a+b{\mathrm{e}}^{\left(-\frac{c}{sin{\alpha }_s}\right)}\right]\hfill \\ a=0.4237-0.00821{\left(6-H\right)}^2\hfill \\ b=0.5055+0.00595{\left(6.5-H\right)}^2\hfill \\ c=0.2711+0.01858{\left(2.5-H\right)}^2\hfill \end{array}\right.$$

(20)

where $G_0$ represents the solar constant, which is valued at 1.366 kW/m², and H denotes the altitude, measured in kilometers.

3. Results and Discussion

To calculate the annual average output thermal power, it is necessary first to determine the average monthly output thermal power. To simplify the calculations, we assume the 21st day of each month represents the average monthly output thermal power. On the 21st of each month, the average output thermal power is determined by taking the mean of the output thermal power at five specific time points: 9:00, 10:30, 12:00, 13:30, and 15:00, which serves as the daily average output thermal power. The output thermal power for each of these moments can be derived from Equation 20. Through this step-by-step approach, the annual average output thermal power can be calculated.

For the cosine efficiency, we take the cosine efficiency for each of the four moments of the heliostat field at the summer solstice, June 22 of that year, and show them as follows:

Figure 4. Summer solstice cosine efficiency

Figure 4. Summer solstice cosine efficiency

From the analysis of the above figure, it can be seen that for different moments of the day, due to the difference between the sun’s altitude angle and the sun’s azimuth angle, it will result in different scope of the mirror field in different moments of the cosine efficiency of the region where the higher cosine efficiency is located, but on the whole, the cosine efficiency of the northern heliostat is larger than that of the southern.

For the collector truncation efficiency for each of the four time points in the heliostat field for the summer solstice, June 22 of that year, the results are shown below:

Figure 5. Summer solstice collector truncation efficiency

Figure 5. Summer solstice collector truncation efficiency

Analyzing the above figure, it can be seen that for the collector truncation efficiency, most of the mirrors in the field of heliosta mirrors are higher, only a small portion of them will show a certain degree of regularity to reduce with the change of the sun’s position, and for this small portion, the closer to the absorption tower, the more obvious the reduction.

For the atmospheric transmittance of each heliostat in this heliostat field, the results are shown below:

Figure 6. Summer solstice atmospheric transmittance efficiency

Figure 6. Summer solstice atmospheric transmittance efficiency

Analyzing the above figure, it can be seen that the atmospheric transmittance in the mirror field shows an extremely obvious regular arrangement, the closer to the absorber tower’s heliostat, the greater its atmospheric transmittance, and vice versa the smaller.

Table 1. This is a two-row, six-column table.

Optical Efficiency	Cosine Efficiency	Shadow Shading Efficiency	Truncation Efficiency	Thermal Power Output (MW)	Hermal Output per Unit Area (kW/m2)
0.5533	0.6836	0.9419	0.9673	33.5254	0.5337

Table 1. This is a two-row, six-column table.

Optical Efficiency	Cosine Efficiency	Shadow Shading Efficiency	Truncation Efficiency	Thermal Power Output (MW)	Hermal Output per Unit Area (kW/m2)
0.5533	0.6836	0.9419	0.9673	33.5254	0.5337

Then, in order to verify the accuracy of the mathematical model of the optical efficiency of heliostat field, the collector diameter is selected and its potential effect on the average optical efficiency is analyzed.

By calculating the truncation efficiency of the floating values respectively, and by observing the sensitivity analysis diagram of the truncation efficiency, we find that the shadow generated by the change of the height and diameter of the collector will not have a great impact on the average optical efficiency, and our model has a certain accuracy.

Figure 7. Truncation efficiency sensitivity analysis

Figure 7. Truncation efficiency sensitivity analysis

4. Conclusions

This study constructed a comprehensive mathematical model for calculating the optical efficiency and power output of a heliostat field. The model is based on the variation of the sun’s rays, and Monte Carlo simulations are used for ray tracing, an innovative method that accurately simulates the trajectory of the rays. At the same time, our calculations are entirely based on fundamental mathematical formulas, ensuring the reliability of the model on an innovative basis. In addition, the reliability of the model is verified by sensitivity analysis.The mathematical model we developed is also very extensible and can be applied to different heliostat fields, and also provides a good basis for future consideration of weather changes and mirror deformation.