基于蒙特卡洛算法的定日镜场的优化设计

Optimization Design of Heliostatic Mirror Field Based on Monte Carlo Algorithm

摘要: 针对定日镜场的优化设计问题，本文首先建立阴影挡光效率模型，结合蒙特卡洛与光线追踪法求解阴影遮挡效率；根据太阳高度角、太阳方位角求解余弦效率、法向直接辐射照度等，得出光学效率与定日镜场的输出热功率。其次，建立以单位面积定日镜输出热功率最大为目标的非线性目标规划模型。在给定相邻定日镜间距的范围内，求得定日镜场中每一圈层的定日镜数目与相邻定日镜的间距大小；采用Campo布置方法确定定日镜坐标，利用阴影挡光效率模型，当额定年平均输出热功率达到60 MW，求得吸收塔位置坐标(0, 0, 80)，定日镜尺寸(宽 × 高)为(6 m × 6 m)，安装高度为4 m，定日镜总数为2475面，定日镜场总面积为89,100 m²。利用Campo方法确定定日镜场中的定日镜排列方式，能够在简化定日镜场的条件下更为简单的求解定日镜的具体坐标。最后，保持额定功率不变，设计定日镜尺寸为(7 m × 7 m)、安装高度为4 m来达到更大单位面积输出热功率。

Abstract: Aiming at the optimization design problem of heliostat field, in this paper, the shadow blocking efficiency model is first established, which combines with Monte Carlo and ray tracing method to solve the shadow blocking efficiency. According to the solar altitude Angle and solar azimuth Angle, cosine efficiency and direct irradiance of normal phase are solved, get the optical efficiency and the output thermal power of heliostatic field. Secondly, a nonlinear target programming model with the maximum output thermal power of heliostat per unit area is established. Within the range given the spacing of adjacent heliostats, the number of heliostats in each circle and the distance between the adjacent heliostats in the heliostats field are obtained. The coordinates of heliostat were determined by Campo’s arrangement method. When the rated annual average thermal output power reached 60 MW, using the shadow blocking efficiency model, the position coordinate of the absorber is (0, 0, 80), and the heliostat size (W × H) is (6 m × 6 m), the installation height is 4 m, the total number of heliostats is 2475, and the total number of heliostats field is 89,100 m². Using Campo method to determine the heliostat arrangement in heliostat field, the specific coordinates of heliostat can be solved more easily under the condition of simplifying heliostat field. Finally, it keeps the rated power unchanged, designs helioscope size of (7 m × 7 m), installation height of 4 m to achieve a large unit area of thermal output power.

1. 引言

塔式太阳能光热发电是我国为了实现“碳达峰”及“碳中和”推出的一种低碳环保的新型清洁能源技术。一方面，系统利用太阳能进行发电，减少二氧化碳等温室气体的排放，减轻空气污染和水污染，保护生态环境。另一方面，作为以新能源为主体的新型电力系统，该技术能够高效利用太阳能、并可以持续进行发电 [1] ，是一种具有巨大潜力的可再生能源技术。本文采用了蒙特卡洛与光线追踪法对定日镜产生阴影面积进行建模，在给定的额定功率下，建立以单位面积定日镜输出热功率最大为目标的非线性目标规划模型。在控制定日镜的各项参数下，完成了对定日镜场的优化设计。本文采用了2023年高教社杯全国大学生数学建模竞赛A题(https://dxs.moe.gov.cn/zx/a/hd_sxjm_sthb/230523/1840580.shtml)进行分析和研究。

2. 模型建立与求解

2.1. 阴影挡光效率模型

定日镜场中太阳光线的损失主要为三个方面：1) 前排定日镜对后排定日镜所接收到太阳光线造成阻挡；2) 后排定日镜反射的太阳光线被前排定日镜遮挡；3) 吸收塔在太阳照射下形成的阴影对部分定日镜造成遮挡。

建立定日镜a镜面坐标系 $\left(X_a,Y_a,Z_a\right)$ 、地面坐标系 $\left(X_d,Y_d,Z_d\right)$ 和定日镜b镜面坐标系 $\left(X_b,Y_b,Z_b\right)$ ，判断定日镜a的反射光线是否落入定日镜b内，确定其在定日镜b的坐标值 $K_2$ ，定日镜间太阳光线路径如图1所示：

Figure 1. Diagram of the path of sun rays between heliostats

图1. 定日镜间太阳光线路径示意图

图中给出某定日镜E中心在地面坐标系的坐标 $\left(X_{e,d},Y_{e,d},Z_{e,d}\right)$ ，定日镜在整个镜场中方位角 ${\theta }_f=\{\begin{array}{l}\arccos \frac{X_{e,d}}{\sqrt{X_{e,d}^2+Y_{e,d}^2}},Y_{e,d}\ge 0\\ \arctan \frac{Y_{e,d}}{X_{e,d}}+\pi ,Y_{e,d}<0\end{array}$ ，太阳入射光线与定日镜b平面法线的夹角 $\theta$ ，太阳高度角 ${\alpha }_s$ [2] 、太阳方位角 ${\gamma }_s$ [3] 的关系为： $\cos \theta =\frac{\sqrt{2}}{2}{\left[\sin {\alpha }_s\cos \beta -\cos \left({\theta }_f-{\gamma }_s\right)\cdot \cos {\alpha }_s\sin \beta +1\right]}^{\frac{1}{2}}$ 。定日镜平面中心法线的俯仰角 $F_y=\arctan \left(\frac{\sin {\alpha }_s+\cos \beta }{2\cos \theta }\right)$ 和方位角 $F_w=\arctan \left(\frac{\sin {\theta }_f\sin \beta -\sin {\gamma }_s\cos {\alpha }_s}{\cos {\theta }_f\sin \beta -\cos {\gamma }_s\cos {\alpha }_s}\right)$ [4] 。

在定日镜a坐标系中的入射点 $K_1=\left(x_a,y_a\right)$ ，经不同坐标系的坐标转换，判断反射光线在定日镜b中是否存在入射点 $K_2$ 。定日镜a坐标系的三个轴在地面坐标系中的向量为 $\left(r_x,s_x,t_x\right),\left(r_y,s_y,t_y\right),\left(r_z,s_z,t_z\right)$ ，转换关系矩阵表达式(2.1) [5] ：

$G=\left[\begin{array}{ccc}{r}_x& r_y& r_z\\ s_x& s_y& s_z\\ t_x& t_y& t_z\end{array}\right]=\left[\begin{array}{ccc}-\sin F_y& -\sin F_w\cos F_y& \cos F_w\cos F_y\\ \cos F_y& -\sin F_w\sin F_y& \cos F_w\sin F_y\\ 0& \cos F_w& \sin F_w\end{array}\right]$ (2.1)

锥形光线在镜面坐标系向量为 $P_j$ ，在地面坐标系向量 $P_0=G\cdot P_j$ 。经过以下5步得到 $K_2=\left(x_b,y_b\right)$ ：

1) 将定日镜a的某一点 $K_1$ 转换为地面坐标系中 ${K^{\prime }}_1$ 。

2) 将地面坐标系的 ${K^{\prime }}_1$ 转换为定日镜b坐标系中的 ${K^{″}}_1$ 。

3) 将地面坐标系的光线向量 $P_0$ 转换为定日镜b坐标系中的 $P_j$ 。

4) 在定日镜b坐标系中计算光线与该镜的交点，由 ${K^{″}}_1=\left({x^{″}}_a,{y^{″}}_a,{z^{″}}_a\right)$ 、 $P_j=\left(x_0,y_0,z_0\right)$ ，解出 $K_2=\left(x_b=\left(z_0{x^{″}}_a-x_0{z^{″}}_a\right)/z_0,y_b=\left(z_0{y^{″}}_a-y_0{z^{″}}_a\right)/z_0,0\right)$ 。

5) 判断 $K_2$ 是否落入定日镜b平面内。

2.2. 大气透射率

定日镜中心到集热器中心距离为 $d_{HR}$ ，大气透射率的计算公式为 [6] ：

${\eta }_{at}=0.99321-0.0001176d_{HR}+1.97\times {10}^{-8}\times d_{HR}^2\left(d_{HR}\le 1000\right)$ (2.2)

2.3. 余弦效率

建立定日镜反射太阳光线图，如图2，集热器中心A，定日镜中心 $O_z$ ，太阳入射光线 $B{O}_z$ ，定日镜反射光线 $O_{z}A$ ，xoy平面的垂线EF和BC，垂足F、C。第i个定日镜入射角 ${\theta }_i$ ，得到余弦效率 [7]

${\eta }_{cos}=\frac{{\left(d_{HR}+\sin {\alpha }_s{h}_a+x_z\cos {\alpha }_s\sin {\gamma }_s+y_z\cos {\alpha }_s\sin {\gamma }_s\right)}^{1/2}}{{\left(2d_{HR}\right)}^{1/2}}$ (2.3)

Figure 2. Three-dimensional representation of sunlight reflected by a heliostat

图2. 定日镜反射太阳光的三维示意图

2.4. 单个定日镜的光学效率

单个定日镜的光学效率 $\eta ={\eta }_{sb}{\eta }_{cos}{\eta }_{at}{\eta }_{trunc}{\eta }_{ref}$ 中阴影遮挡效率 ${\eta }_{sb}$ 、余弦效率 ${\eta }_{cos}$ 和大气透射率 ${\eta }_{at}$ 由之前式子已表示出，资料显示集热器截断效率 ${\eta }_{trunc}$ 取值60%~90% [8] ，镜面反射率 ${\eta }_{ref}$ 为0.92。

2.5. 法向直接辐射辐照度

法向直接辐射辐照度 [9] 为

$DNI=G_0\left[a+b\exp \left(-\frac{c}{\sin {\alpha }_s}\right)\right]$

其中： $a=0.4237-0.008216{\left(6-H\right)}^2$ ， $b=0.5055+0.00595{\left(6.5-H\right)}^2$ ， $c=0.2711+0.01858{\left(2.5-H\right)}^2$ ， $G_0=1.366\text{kW}/{\text{m}}^{\text{2}}$ 是太阳常数，H为当地的海拔高度为3000 m。

2.6. 定日镜场年平均输出热功率的求解

每面定日镜采光面积 $A_i$ 即定日镜尺寸为6 m × 6 m，由之前给出的各个定日镜的光学效率 ${\eta }_i$ 和定日镜场的法向直接辐射辐照度DNI，定日镜总数N为1745面，按照以下公式即可求出定日镜场的输出热功

率： $E_{field}=DNI\cdot {\displaystyle \underset{i}{\overset{N}{\sum }}{A}_i{\eta }_i}$ ，得出平均光学效率 $\bar{\eta }$ ，平均余弦效率 $\overline{{\eta }_{cos}}$ ，平均阴影遮挡效率 $\overline{{\eta }_{sb}}$ ，平均截断效率 $\overline{{\eta }_{trumc}}$ ，单位面积镜面平均输出热功率 $\overline{\left(\frac{E_{field}}{N\times A}\right)}$ ，结果如表1和表2所示：

Table 1. Average optical efficiency and output power on 21st per month

表1. 每月21日平均光学效率及输出功率

Table 2. Annual average optical efficiency and output power

表2. 年平均光学效率及输出功率

2.7. Campo布置方法

以吸收塔为中心构造定日镜场为同心圆模型，每个定日镜的占地范围视为圆形，依据Campo布置方法布置圆形定日镜场，以定日镜间不发生机械碰撞为原则进行周向均匀布置，模型示意图如图3：

Figure 3. Heliostat field model

图3. 定日镜场模型

假定定日镜位置坐标为 $\left(X,Y\right)$ ，其中 $\left(x_i,y_i\right)$ 表示定日镜场中第i个定日镜在地面坐标系下的坐标。 $\left(x_{b,i},y_{b,i}\right)$ 表示光线在定日镜b中的反射点坐标， $x_{i,j}$ 表示第j圈圆周上的第i个点。构建非线性规划模型表达式如下：

$\begin{array}{l}{E}_{field}=DNI{\displaystyle \underset{i=1}{\overset{N}{\sum }}{A}_i{\eta }_i}\\ \{\begin{array}{l}DNI=G_0+\left[a+b\exp \left(-c/\sin {\alpha }_s\right)\right]\\ {\eta }_i={\eta }_{sb,i}{\eta }_{cos,i}{\eta }_{at,i}{\eta }_{trunc,i}{\eta }_{ref,i}\\ {\eta }_{at,i}=0.99321-0.0001176d_{HR}+1.79\times {10}^{-8}{d}_{HR}^2\left(d_{HR}\le 1000\right)\\ {\eta }_{cos,i}=\cos {\theta }_i=d_{HR}+h_a\sin {\alpha }_s+x_i\cos {\alpha }_s\sin {\gamma }_s+y_i\cos {\alpha }_s\sin {\gamma }_s/{\left(2d_{HR}\right)}^{\frac{1}{2}}\\ {\eta }_{sb,i}=\left(x_{b,i}-3\right)\left(y_{b,i}-3\right)/36\\ 100<\left|x_i\right|<350\\ 100<\left|y_i\right|<350\\ 11<a<14\\ \sqrt{{\left(x_{i+1,j}-x_{i,j}\right)}^2+{\left(y_{i+1,j}-y_{i,j}\right)}^2}>a\\ \sqrt{{\left(x_{i,j}-x_{i,j+1}\right)}^2+{\left(y_{i,j}-y_{i,j+1}\right)}^2}>a\\ r_j=\{\begin{array}{l}\begin{array}{cc}106& j=1\end{array}\\ \begin{array}{cc}106+aj& j>1\end{array}\end{array}\left(j表示的是第j圈;其中100<\left|r_j\right|<350\right)\end{array}\end{array}$ (2.4)

利用Python进行编程，将定日镜间的距离从11 m开始，以0.1 m的步长增加至14 m，R表示定日镜场的圆周半径，r为定日镜场圆周上单个定日镜中心的半径，依据定日镜之间距离以及R求得每个圈层内设置的定日镜个数 $n=360\times \left[2\times \arcsin \left(r/R\right)\right]$ 。

编程所绘制的不同间距定日镜安排示意图4如下：

Figure 4. Schematic diagram of heliostat arrangement with different spacing

图4. 不同间距定日镜安排示意图

得出最优的相邻定日镜之间的间距为12 m，在满足定日镜场的额定年平均输出热功率为60 MW的情况下，求得相关数据为：吸收塔位置坐标(0, 0, 80)，定日镜尺寸(宽 × 高)为(6 m × 6 m)，其安装高度为4 m，定日镜总数为2475面，定日镜场总面积为89,100 m²。

2.8. 求解更大的单位面积输出热功率

保持额定功率不变，设计定日镜尺寸和安装高度达到更大的单位面积输出热功率。现将定日镜的尺寸从原来的6 m × 6 m增加至7 m × 7 m [10] ，使相邻镜面在不超过间距范围的同时达到最大的镜面尺寸；接着为了尽可能增大定日镜的输出热功率，使定日镜从内圈至外圈的安装高度逐次递增，由题可知，安装高度的范围是3.5 m~6 m。

为确定各圈层所包含的定日镜数目，利用MATLAB进行编程 [11] ，搜索满足各圈层半径范围内定日镜的个数，再逐个相减，即可求得从内至外各圈层内包含的定日镜的数目，之后将安装高度分成20等份，在Excel表格中为各圈层内的定日镜附加上其安装高度的数值。最后将所得的数值代入上述的Python代码，并增加上关于定日镜安装高度的关系函数，即可求得各项数值。结果如表3、表4和表5所示：

Table 3. Average optical efficiency and output power on the 21st day of each month after optimization

表3. 优化后的每月21日平均光学效率及输出功率

Table 4. Annual average optical efficiency and output power

表4. 年平均光学效率及输出功率

Table 5. Design parameter

表5. 设计参数

3. 结论

本文针对定日镜场的优化设计问题，建立阴影挡光效率模型和以单位面积定日镜输出热功率最大为目标的非线性目标规划模型。

在定日镜尺寸确定情况下，得出表1和表2定日镜场的平均年光学效率、年平均输出热功率，以及单位镜面面积年平均输出热功率。

在定日镜场额定年平均输出热功率为60 MW情况下，得出最优的相邻定日镜之间的间距为12 m，在满足定日镜场的额定年平均输出热功率为60 W的情况下，求得相关数据为：吸收塔位置坐标(0, 0, 80)，定日镜尺寸(宽 × 高)为(6 m × 6 m)，其安装高度为4 m，定日镜总数为2475面，定日镜场总面积为89,100 m²。

在定日镜场达到额定功率的条件下，单位镜面面积年平均输出热功率尽量大，得出按表5的设计参数。

NOTES

^*共同第一作者。