A Parameter Optimization Design and Numerical Calculation Method for Circular Heliostat Field

Solar energy is a low-carbon, eco-friendly, clean energy source that is easily converted into heat and electricity, making the study of solar thermal power generation technology of great significance and with broad application prospects. To improve the performance of solar power plants, it is necessary to reasonably arrange the light reflectors. Both denser heliostats are needed to ensure enough light can be concentrated into the solar collector, thereby improving the power of the power station; at the same time, waste of material must be prevented, striving to increase output thermal power per unit mirror area. This paper is based on the background of the 2023 CUMCM Problem A. Firstly, it analyzes parameters such as the optical efficiency of the original heliostat field, the annual average output power, and the annual average thermal power output per unit mirror area. Secondly, based on the data provided in the appendix, it analyzes the distribution pattern of the heliostat field and establishes an optimization design model for the heliostat field. Next, the optimization design model of the heliostat field is solved using discrete numerical calculation methods. Finally, according to the results of the model solution, the placement positions of the heliostats, the size of the mirrors, and the installation height are adjusted to improve the annual average thermal power output per unit mirror area, thereby achieving an optimized design of the heliostat field.


Introduction
The plan is to construct a tower-type solar thermal power station [1], he heliostat field is located at 98.50 degrees east longitude and 39.40 degrees north latitude, with an altitude of 3000 meters, it is arranged within a circular area with a radius of 350 meters.The center of the circular area is selected as the origin, the positive direction of the east is taken as the positive direction of the x-axis, the positive direction of the north as the positive direction of the y-axis, and the upward direction perpendicular to the ground as the positive direction of the z-axis to establish the mirror field coordinate system.The absorption tower at the center of the heliostat field is 80 meters high, the diameter of the cylindrical external light-receiving heat collector is 7 meters and its height is 8 meters.No heliostats are installed within 100 meters of the absorption tower, the space is reserved for building factories to install power generation, energy storage, and control equipment.Each piece of heliostat is 6 meters in height and width, installed at a height of 4 meters, and the planned installation location of each heliostat is known [1].
This article first calculates the optical efficiency, average annual output power, and average annual thermal output power per mirror area in existing tower solar thermal power stations; then, by adjusting radians per hour.Assuming that the Earth rotates uniformly if the local time is , then the current solar hour angle  is a linear function of  [2].

Sun Declination Angle 𝛿𝛿
Sun declination is the angle between the plane of the Earth's equator and the line connecting the center of the sun and the earth.That is, the latitude value of the solar zenith, which moves within the range of +23°26′ and -23 °26′, becomes a symbol of the seasons.With March 21st, the vernal equinox, as the 0th day, and  set as the accumulated day, the calculation formula for the solar declination angle is [3]:

Solar altitude Angle 𝛼𝛼 𝑠𝑠
The angle between the line  connecting the Sun  and the origin, and the  plane, is called the solar altitude angle α s , where 0 ≤   ≤  2 . If the local latitude is , then the calculation formula for the solar altitude angle is [4]: sin  = cos  cos  cos  + sin  sin . (3)

Sun Azimuth Angle 𝛾𝛾 𝑠𝑠
In the mirror field coordinate system, the angle between the vector  �����⃗ and the positive direction of the -axis or true north is known as the sun azimuth angle   , − ≤   ≤ , and it is conventionally positive for clockwise rotation.The calculation formula is [5]: (4)

The Normal Direct Radiation Irradiance of Solar Rays 𝐷𝐷𝐷𝐷𝐷𝐷
The solar radiation energy received per unit area per unit time on a plane perpendicular to the sunlight on the earth is called the normal direct radiation irradiance, denoted as , it is related to the solar elevation angle and the local altitude  (unit in kilometers), which can be calculated by the following approximate formula [6] thereinto,  = 0.4237 − 0.00821(6 − ) 2 ,  = 0.5055 + 0.00595(6.5− ) 2 ,  = 0.2711 + 0.01858(2.5− ) 2 ,  0 = 1.366/ 2 . (6) 2.6.Cosine Efficiency   Incident light  and the normal direction n of the mirror reflection point must have an angle , which will cause the light to be obliquely spread on the mirror surface, thereby converting the energy of the light from being distributed on the head-on area to being obliquely spread on the mirror surface, which will cause power loss per unit area, called cosine loss; after deducting the cosine loss, the ratio of the remaining energy to the original light energy is called the cosine efficiency, denoted as Let the average Earth-Sun distance be   , the Earth's radius be   , the altitude of the heliostat field be , and the latitude be .Let the accumulated day be  (in days), and the local time be  (in hours).In the coordinate system of the mirror field, set the center coordinates of a certain heliostat as (, , ).The height of the absorption tower is marked as   , so the coordinates of the collector center are (0,0,   ).
Connect the origin  and the point light source  of the sun to get the vector  �����⃗ , the modulus of which is: � �����⃗ �=R  −   −  , according to the conversion between spherical coordinates and Cartesian coordinates, the Cartesian coordinates of the sun in the mirror field coordinate system are S(  ,   ,   ), thereinto, (9)

Atmospheric Transmissivity 𝜂𝜂 𝑡𝑡𝑡𝑡
According to literature [7], during the process of light collection and reflection in a heliostat, the light is scattered by the atmosphere or obstructed by particles in the air, resulting in energy loss.The remaining light energy after loss, as a proportion of the energy before loss, is called the atmospheric transmittance   , it is related to factors such as the reflection distance   (in meters), atmospheric pressure, humidity, air purity, etc.This paper uses an approximate formula [8]   = 0.99321 − 0.0001176  + 1.97 × 10 −8 ×   2 (  ≤ 1000). (10)

Collector Cutoff Efficiency 𝜂𝜂 𝑡𝑡𝑡𝑡𝑡𝑡𝑛𝑛𝑐𝑐
Because the shape of the collector is not infinitely large, only a part of the light reflected by the heliostat is received by the collector, and some of the light leaks out from the side of the collector.Therefore, the collector only receives reflected light that meets certain geometric conditions.The truncation efficiency   of the collector can be defined as the ratio between the received light and the total reflected light According to the law of plane mirror imaging, the light reflected by the heliostat, when extended backwards, must intersect at point  ′ , the symmetrical point of the sun in the plane mirror.Because vectors  ������⃗ and  ′  �������⃗ are symmetrical about the heliostat, we have � ′  �������⃗ � = � ������⃗ �, which determines the magnitude of  ′  �������⃗ .Then, because points ′ ,  and  J are collinear, we can determine the direction of vector  ′  �������⃗ .This gives us the result Furthermore, we obtain the coordinates of vector  ′ �������⃗ , which also equals to the coordinates of point  ′ : Next, consider the projection of the rectangular sundial illuminated by the point light source  ′ on the  plane.The top and bottom of the sundial always remain horizontal, so the direction vector of the top can be set as {1, , 0}.Then, based on the top being perpendicular to  �⃗, it is determined that  = −     , set the upper edge length of the sundial as ℎ and ℎℎ, and derive the upper edge vector as: The direction vector of the left side of the sundial can be obtained from the vector product of the top vector and the normal vector: So the left vector is: Then, calculate the length of the positive projection on the reflected ray  �����⃗ of the upper frame of the sundial.
Similarly, the length of the positive projection on the reflected ray  �����⃗ of the lower frame of the sundial is Similarly, consider the image that  ′ projects the collector onto the  plane is the positive projection of the rectangular cross-section of the cylindrical collector along the reflected light ray  �����⃗ .Due to symmetry, the width of the cross-sectional rectangle is the diameter of the collector, 7 meters, and the height of the rectangle is the length of the positive projection of the collector height vector  = {0,0,8} on the reflected light ray  �����⃗ : By symmetry, it is known that the centers of the two projected rectangles coincide, and the effective receiving area of the collector is the overlap area of the two rectangles.As shown in the figure 2，the width of the overlapped rectangle is min{ ′ , 7}, the height is min{ ′ ,  ′ }, and its area is min min{ ′ , 7} ⋅ min{ ′ ,  ′ } Hence, the truncation efficiency of the collector is as follows: 2.9.Mirror Reflectivity   It's inevitable that sundial mirrors will accumulate dust or impurities, and they are also subject to erosion from rain, snow, and frost.The ratio of the energy of reflected light to the energy of incident light is known as the mirror reflectivity   , which generally ranges between 0.93-0.94.Here it can be taken as 0.935.

Optical Efficiency of Heliostat 𝜂𝜂
For the sake of simplicity, assume that various factors causing light loss, such as mirror reflectivity, collector truncation efficiency, and atmospheric transmissivity, are independent of each other.Under these independence assumptions, each influencing factor does not interfere with each other, so we can discuss their affects on light loss separately.The optical efficiency of the sundial mirror can then be obtained by the probability multiplication formula,  =  cos ⋅   ⋅   ⋅   .

Thermal Output Power of Heliostat Field
The solar radiant energy received on a unit area perpendicular to the sunlight on the earth within a unit time is ; for the  heliostat, denote its light-collection area as   , then the solar radiant energy it receives is  ⋅   ; Then consider the cosine efficiency, atmospheric transmittance, collector truncation efficiency, mirror reflectance and a series of factors causing light energy loss resulting in the total factor   ,The contribution of the  heliostat to the thermal power of the collector is:  ⋅   ⋅   .By traversing N heliostats and summing them up, we get the output thermal power of the heliostat field at the current moment: . (20)

The Average Annual Output Thermal Power of Heliostat Field
Theoretically, to calculate the annual average output heat power of a heliostat field, one should use the infinitesimal method to divide the year into numerous time elements [,  + ] .Within this infinitesimal time period, the sunlight doesn't have time to change.
Thereby the heliostats within the field remain stationary, allowing us to calculate the output thermal energy element of the heliostat field.
Then, integrating with respect to , we obtain the annual output thermal energy of the heliostat field  = �   () ⋅  ∈ .Hence, we can then derive the annual average output heat power of the To avoid complex calculus operations, a numerical integration method can be adopted.It only requires obtaining the instantaneous output thermal power of the heliostat field at several sampling points throughout the year.Then, by taking the weighted average of these function values at the sampling points, we can get an approximation of the integral, which is the annual average output thermal power of the heliostat field.
To further simplify, this paper only calculates the annual average output thermal power of the heliostat field at 60 sampling points: on the 21st of each month throughout the year (12 months in total) and five representative times each day (9 a.m, 10:30 a.m, 12 p.m, 1:30 p.m, 3 p.m).The arithmetic mean of these 60-time sampling points is taken as the crude approximation of the annual average output thermal power of the heliostat field: Furthermore, by dividing the annual average output thermal power of the heliostat field by the sum of the mirror area of the heliostat, we get the annual average output thermal power per unit mirror area.

Optical Efficiency Analysis of the Original Planned Heliostat Field
By substituting the relevant parameters of the original planned heliostat field, we set mirror reflectivity   = 0.935, average distance between the Sun and Earth as   = 149590000 kilometers, Earth's radius   = 6371.393kilometers, the altitude of the heliostat field = 3000 meters, latitude  = 39.4°,height of the absorption tower   = 80 meters, and ℎ = 6 and ℎℎ = 6 of the heliostat.We then calculate the efficiency and other physical quantities of the heliostat field in sequence.The results are shown in the following table:   As shown in the figure 3, it can be seen that the heliostats are distributed on a series of concentric circles.Secondly, we transformed the center coordinates of the heliostat into polar coordinates and extracted the polar radius array.It was found that there were 18 tracks in total, numbered from inside to outside.We plotted the relationship diagram of track radius with respect to track number, which is very close to a straight line.The linear correlation coefficient is 0.999999999998667, indicating that the radius can be considered as an arithmetic sequence.According to the linear fitting equation, the common difference of the track radius is about 13.5.
By extracting the polar angles of the center points of the heliostats on each layer of track, it was found that the polar angles on the same track are approximately arranged in arithmetic progression.For tracks with odd numbers, the first term of their polar angle is zero; for tracks with even numbers, the first term of their polar angle is about half of the common difference of the polar angle.The common difference of the polar angles on each track decreases with the increase of track number, and its common difference is approximately -0.003.
Therefore, we can summarize the design scheme of the heliostat field as: Select the first track 100 meters away from the origin, the radius of which is recorded as  1 (for example, the radius of the first track selected in [1] is 107).The polar angles are arranged in arithmetic progression, with the first term being 0, the angular tolerance is recorded as  1 , and a total of  1 = � 2  1 � heliostats are arranged; Let the tolerance of the track radius be   (the tolerance in [1] is about 13.5), and the angular tolerance on the second track is smaller than the angular tolerance on the first track by  ( = 0.003 in [1]), then the radius of the second track is  2 =  1 +   , and the angular tolerance on the second track is  2 =  1 − , a total of  2 = � 2  2 � heliostats are arranged, the first term of the polar angle is

Figure 2 .
Figure 2. Schematic diagram of effective receiving area of the collector.

4. 1 .
The Field Distribution Law of Heliostat For proper design of the heliostat field, we first create the projection diagram of the heliostat center on the ground plane [1].

Table 1 .
Annual average optical efficiency and output power of the original planned heliostat field.