定日镜场光学效率建模:从物理原理到Python实现
2026/8/27 4:13:36 网站建设 项目流程

1. 从赛题到模型:一次完整的建模实战复盘

去年带队参加数学建模竞赛,我们组选的就是这道“定日镜场的优化设计”。说实话,当时看到题目,尤其是问题一,感觉既兴奋又棘手。兴奋在于,这是一个典型的工程优化问题,有明确的物理背景和现实意义;棘手在于,它不像一些纯数据分析题那样有现成的数据集,一切都要从物理原理和几何关系出发,自己“无中生有”地构建模型。问题一的核心,说白了,就是给定一个特定的时刻(比如某年某月某日的中午12点),在已知太阳位置、接收塔位置和高度、以及单个定日镜尺寸与安装高度的情况下,如何计算镜场中任意一面定日镜的“光学效率”。这个效率,直接决定了这面镜子在那一刻能往塔顶接收器上反射多少有效太阳能,是整个镜场优化设计的基石。如果你也在准备类似的竞赛,或者对太阳能光热发电的建模感兴趣,那么跟着我一起拆解这个问题,你会看到如何将一道看似抽象的赛题,转化为一步步可计算、可编程的数学模型。这个过程,远比直接套用一个现成的公式更有价值。

2. 问题一的核心:拆解“光学效率”的计算链条

拿到问题一,首要任务不是急着写公式,而是彻底理解“光学效率”这个目标量到底由哪些因素决定。根据题目描述和光热发电的基础知识,一面定日镜的光学效率(η_optical)通常不是单一值,而是几个子效率的乘积。我们的模型建立,本质上就是为每一个子效率找到准确的数学表达。经过文献调研和小组讨论,我们将其分解为以下四个主要部分,这也构成了我们模型的核心框架:

2.1 余弦效率 (η_cos):入射角带来的能量损失

这是最直观的一个效率因子。定日镜的法线方向需要精确调整,使得入射的太阳光经反射后正好指向塔顶的接收器。然而,太阳光线并非垂直照射镜面,而是存在一个入射角。根据兰伯特余弦定律,镜面实际接收到的太阳辐射强度与入射角余弦(cosθ)成正比。当太阳光垂直入射时(θ=0°),cosθ=1,效率最高;入射角越大,有效接收面积越小,效率越低。 因此,余弦效率 η_cos = cos(θ_i),其中 θ_i 是太阳光线与定日镜法线之间的夹角。计算这个夹角,就需要知道三个关键向量的方向:太阳光线方向向量、定日镜中心点到塔顶接收器中心的向量(反射光方向)、以及定日镜的法线方向。根据光学反射定律(入射角等于反射角,且入射光线、法线、反射光线共面),法线方向恰好是太阳光线方向向量和反射光线方向向量的角平分线方向。由此,我们可以通过向量运算精确求出 θ_i。

注意:这里容易混淆“入射角”和“太阳高度角/方位角”。入射角是相对于镜面局部坐标系的概念,必须通过空间向量计算,不能直接用太阳的位置角替代。

2.2 阴影遮挡效率 (η_shadow):镜子间的“互相伤害”

在一个密集布置的镜场中,镜子之间难免会互相遮挡阳光,或者遮挡反射光路。这主要分为两种:

  1. 阴影损失:位于前排的定日镜挡住了太阳光,使其无法照射到后排镜子的部分或全部镜面。
  2. 遮挡损失:位于前排的定日镜挡住了后排镜子反射向接收塔的光路。 在问题一中,通常假设镜场规模不大或镜子间距合理,且计算的是特定瞬时效率,有时可以忽略阴影遮挡,即设 η_shadow = 1。但如果题目要求考虑,或者作为模型完备性的一部分,就需要进行几何判断。我们需要计算太阳光线方向下,镜面A的投影是否会覆盖镜面B;以及从镜面B中心到塔顶的连线上,是否会与镜面A相交。这涉及到三维空间中的多边形投影和光线相交检测,计算量较大,是编程实现中的一个难点。

2.3 大气透射效率 (η_atm):光在空气中走过的“损耗”

太阳光从定日镜反射到塔顶接收器的过程中,需要穿越一段大气距离。大气中的尘埃、水汽等会吸收和散射部分光能,导致能量衰减。这种衰减通常用布格-朗伯定律(Bouguer-Lambert Law)的简化形式来描述:η_atm = exp(-k * D)。其中,D 是定日镜中心到塔顶接收器的空间直线距离(斜距),k 是大气衰减系数,是一个经验常数,题目一般会给定(例如,k=0.0001 m⁻¹ 量级)。这意味着,镜子离塔越远,光路越长,大气衰减越严重。这个计算相对简单,关键在于准确计算每一面镜子到塔顶的三维距离。

2.4 截断效率 (η_trunc):接收器没接住的“溢出”光能

即使反射光路准确指向接收器,由于太阳本身不是一个点光源,而是具有约0.5°张角的盘面,以及镜面并非理想光学平面存在一定的斜率误差,导致反射光斑不是一个理想的点,而是一个有一定大小的光斑。当这个光斑大于接收器的开口尺寸时,就会有一部分光能没有被接收器捕获,而是“溢出”了,这部分损失就是截断损失。 计算截断效率需要建立光斑模型。一个常用的简化方法是“圆锥面”模型:假设反射光束是一个以理想反射光线为轴、具有一定锥角(由太阳形状张角和镜面光学误差共同决定)的圆锥。接收器则被视为一个垂直于地面(或有一定倾角)的平面矩形或圆形区域。截断效率就是该圆锥光束与接收器平面相交部分的光通量占总光通量的比例。这通常需要通过数值积分或蒙特卡洛光线追迹来精确计算,在竞赛限时条件下,可能会采用基于误差分布函数的解析近似公式。 对于问题一,如果题目未强调或提供误差参数,有时会先假设为理想情况,即光斑完全被接收器接收,η_trunc = 1。但更严谨的做法是,即使简化,也应说明这一项的存在和可能的处理方式。

3. 模型建立的关键步骤:从物理到数学公式

理解了效率的构成,下一步就是用数学语言精确描述每一个环节。这里我分享一下我们当时建立的模型框架,你可以把它看作一个可执行的“计算清单”。

3.1 第一步:定义坐标系与关键参数

建立一个清晰的空间坐标系是所有计算的基础。我们采用如下右手坐标系:

  • 原点O:定日镜场所在水平地面上的某一点(例如,场地的西南角或中心)。
  • X轴:指向正东。
  • Y轴:指向正北。
  • Z轴:垂直地面向上。关键参数(题目给定或假设):
  • 接收塔位置:T = (x_t, y_t, h_t),其中h_t为塔高。
  • 定日镜位置:H_i = (x_i, y_i, h_m),其中h_m为定日镜的安装高度(镜面中心离地高度),通常所有镜子相同。
  • 定日镜尺寸:假设为矩形,长L,宽W。
  • 太阳位置:由太阳高度角α_s和太阳方位角γ_s(从正北顺时针起算)确定。这两个角可以根据比赛题目给出的具体日期、时间和地点,通过太阳位置算法(如SPA算法)精确计算,这是另一个子模块。
  • 大气衰减系数:k。
  • 太阳形状张角:δ(约0.0093 rad)。
  • 镜面光学误差(斜率误差):σ_slope(题目可能给定)。

3.2 第二步:计算太阳与反射光线的方向向量

  • 太阳光线单位向量 (S):由太阳高度角和方位角计算。 S = [S_x, S_y, S_z] = [cos(α_s) * sin(γ_s), cos(α_s) * cos(γ_s), sin(α_s)] 注意:这里方位角γ_s的定义需与坐标系一致。我们采用从北顺时针,所以向量分量如此。
  • 反射光线单位向量 (R):从定日镜中心H指向塔顶T。 R = (T - H) / ||T - H||,其中“|| ||”表示向量的模(长度)。

3.3 第三步:计算法线向量与余弦效率

根据反射定律,法线向量 N 是入射光线反向向量(-S)和反射光线向量(R)的角平分线方向。 N = (R - S) / ||R - S|| (注意:这里用R - S,因为入射方向为-S,反射方向为R,它们的和向量方向即角平分线方向,需归一化)那么,入射角 θ_i 即为向量(-S)与N的夹角,或者向量R与N的夹角。cos(θ_i) = |(-S) · N| = |R · N| (因为根据反射定律,这两个点积的绝对值相等) 因此,余弦效率 η_cos = |R · N|。 这里取绝对值是因为我们只关心夹角大小,方向不影响余弦值。同时,这也保证了效率值为正。

3.4 第四步:计算大气透射效率

首先计算定日镜到塔顶的直线距离 D = ||T - H||。 然后,大气透射效率 η_atm = exp(-k * D)

3.5 第五步:截断效率的简化建模

在竞赛有限时间内,实现完整的光线追迹不现实。我们采用了一种基于高斯误差假设的解析近似方法,这也是许多工程简化模型的做法。 假设反射光斑在接收器平面上的能量分布是一个二维圆对称高斯分布。接收器为半径为R_ap的圆形开口。那么,截断效率可以近似为: η_trunc ≈ erf( sqrt(2) * R_ap / (σ_total * D) )^2 其中,erf是误差函数。σ_total 是总的角误差标准差,由太阳张角δ和镜面光学误差σ_slope合成:σ_total = sqrt( (δ/2)^2 + (2σ_slope)^2 )。 这个公式的物理意义是:光斑的扩展角度标准差是σ_total,经过距离D后,在接收器平面上的光斑半径标准差约为σ_total * D。接收器能接收到的光能比例,就是这个高斯分布落在半径为R_ap的圆内的概率。 如果题目未给出误差参数,可暂时设η_trunc=1,但必须在模型说明中阐述此项。

3.6 第六步:综合光学效率模型

最终,单面定日镜在特定时刻的光学效率为:η_optical = η_cos * η_shadow * η_atm * η_trunc对于问题一,在忽略阴影遮挡的情况下,模型简化为:η_optical = η_cos * η_atm * η_trunc至此,我们得到了一个输入(镜子坐标、太阳位置、系统参数)到输出(光学效率)的完整数学模型。这个模型是高度参数化的,改变任何一个输入参数,都能重新计算效率。

4. 模型求解与编程实现:将公式转化为代码

模型建立后,求解就是编程计算。我们当时使用Python,因其科学计算库强大。以下是核心实现步骤和代码片段思路。

4.1 环境准备与太阳位置计算

首先需要能计算太阳位置。我们使用了pysolar库(需注意时区、地理位置输入)作为备用,但比赛中更常用的是自己实现一个简化版的SPA(Solar Position Algorithm)函数,以确保可控性和无依赖。这里假设我们已经有了一个函数get_sun_position(lat, lon, year, month, day, hour, minute, second),它返回太阳高度角(alpha_s)和方位角(gamma_s)。

import numpy as np import math # 系统固定参数(示例值,实际由题目给出) k = 0.0001 # 大气衰减系数 (m^-1) h_t = 80 # 塔高 (m) h_m = 4 # 定日镜安装高度 (m) R_ap = 1.5 # 接收器半径 (m) delta = 0.0093 # 太阳半张角 (rad) sigma_slope = 0.001 # 镜面斜率误差 (rad),示例 # 接收塔位置 (以镜场中心为原点示例) T = np.array([0, 0, h_t]) # 假设有一面定日镜位于 (x_i, y_i) H = np.array([50, 30, h_m]) # 示例坐标 # 假设的日期时间地点(用于计算太阳位置) lat, lon = 40.0, 110.0 # 纬度,经度 year, month, day, hour, minute, second = 2023, 6, 21, 12, 0, 0 # 夏至日正午 # 计算太阳高度角和方位角 (这里调用函数,实际需实现) alpha_s, gamma_s = get_sun_position(lat, lon, year, month, day, hour, minute, second) alpha_s_rad, gamma_s_rad = np.radians(alpha_s), np.radians(gamma_s)

4.2 核心效率计算函数实现

接下来,实现3.2到3.5节的各个效率计算函数。

def calculate_cosine_efficiency(sun_vec, target_vec): """ 计算余弦效率 sun_vec: 太阳光线单位向量 (指向太阳) target_vec: 定日镜到目标的单位向量 (指向接收塔) """ # 计算法线向量 (角平分线方向) normal_vec = (target_vec - sun_vec) normal_vec = normal_vec / np.linalg.norm(normal_vec) # 余弦效率 = |反射光线·法线| = |目标向量·法线| eta_cos = abs(np.dot(target_vec, normal_vec)) # 理论上入射角不应超过90度,这里可加个约束 eta_cos = max(0, min(1, eta_cos)) return eta_cos, normal_vec def calculate_atmospheric_efficiency(distance, k_coeff): """计算大气透射效率""" return math.exp(-k_coeff * distance) def calculate_truncation_efficiency(distance, R_ap, delta, sigma_slope): """计算截断效率 (高斯近似)""" # 总角误差标准差 sigma_total = math.sqrt((delta/2)**2 + (2*sigma_slope)**2) # 参数 u = R_ap / (sigma_total * distance), 避免除零 if distance == 0 or sigma_total == 0: return 1.0 u = R_ap / (sigma_total * distance) # 误差函数近似计算 (可使用math.erf,或数值近似) from math import erf eta_trunc = erf( math.sqrt(2) * u )**2 return eta_trunc def calculate_optical_efficiency(H, T, sun_alpha, sun_gamma, k, h_m, R_ap, delta, sigma_slope): """ 计算单面定日镜的光学效率 H: 定日镜中心坐标 [x, y, z] T: 接收塔顶坐标 [x, y, z] sun_alpha: 太阳高度角 (弧度) sun_gamma: 太阳方位角 (弧度,从北顺时针) """ # 1. 计算太阳方向向量 (指向太阳) S = np.array([ math.cos(sun_alpha) * math.sin(sun_gamma), math.cos(sun_alpha) * math.cos(sun_gamma), math.sin(sun_alpha) ]) # 2. 计算定日镜到塔顶的向量和距离 vec_HT = T - H distance = np.linalg.norm(vec_HT) R = vec_HT / distance # 反射光线单位向量 # 3. 计算各项效率 eta_cos, _ = calculate_cosine_efficiency(S, R) eta_atm = calculate_atmospheric_efficiency(distance, k) eta_trunc = calculate_truncation_efficiency(distance, R_ap, delta, sigma_slope) # 4. 综合光学效率 (暂不考虑阴影遮挡) eta_optical = eta_cos * eta_atm * eta_trunc return eta_optical, eta_cos, eta_atm, eta_trunc, distance

4.3 对镜场批量计算与结果分析

对于问题一,通常需要计算镜场中所有镜子在特定时刻的效率。这只需将上述函数放入循环即可。

# 假设 mirrors 是一个列表,包含所有定日镜的 [x, y] 坐标 mirrors = [[50, 30], [60, 20], [40, 40], ...] # 示例数据 h_m = 4 # 镜高 results = [] for x, y in mirrors: H_i = np.array([x, y, h_m]) eta_opt, eta_c, eta_a, eta_t, d = calculate_optical_efficiency( H_i, T, alpha_s_rad, gamma_s_rad, k, h_m, R_ap, delta, sigma_slope ) results.append({ 'position': (x, y), 'eta_optical': eta_opt, 'eta_cos': eta_c, 'eta_atm': eta_a, 'eta_trunc': eta_t, 'distance': d }) # 分析结果,例如找出效率最高/最低的镜子,计算平均效率等 efficiencies = [r['eta_optical'] for r in results] avg_eta = np.mean(efficiencies) max_eta = max(efficiencies) min_eta = min(efficiencies) print(f"镜场平均光学效率: {avg_eta:.4f}") print(f"最高效率: {max_eta:.4f}, 最低效率: {min_eta:.4f}")

通过这样的批量计算,我们就能得到镜场在指定时刻的瞬时性能快照。这为后续问题(如年总输出优化、镜场布局优化)提供了至关重要的基础数据。

5. 建模过程中的关键陷阱与应对策略

在实现上述模型的过程中,我们踩过不少坑,也总结出一些让模型更稳健、计算结果更可靠的经验。

5.1 向量计算中的方向与符号陷阱

这是最容易出错的地方。太阳方向向量S,是定义为“从地面点指向太阳”还是“从太阳指向地面点”?反射向量R,是“从镜子指向塔”还是“从塔指向镜子”?法线向量N的计算公式是**(R - S)还是(S - R)**?

  • 我们的经验:必须严格统一物理定义。我们定义:
    • S:太阳光线方向单位向量,即光传播的方向(从太阳到镜子)。
    • R:反射光线方向单位向量,即光传播的方向(从镜子到塔)。
    • 那么,根据反射定律,镜面法线N应是入射方向(-S)和反射方向(R)的角平分线,即N = normalize(R - S)。因为R - S = R + (-S),正是两个方向向量的和向量,指向角平分线。
  • 验证方法:用一个特例验证。例如,假设太阳在正东,镜子在原点正东10米,塔在原点正上方。此时,S大概为[-1,0,0](假设东为X轴正),R为[0,0,1](指向正上),计算出的N应该大致在X-Z平面的45度方向。用程序算出后,检查点积dot(N, R)dot(N, -S)是否近似相等(即入射角等于反射角余弦值)。

5.2 角度制与弧度制的混乱

三角函数(sin,cos,arctan等)在绝大多数编程语言(包括Python的math/numpy)中默认使用弧度制。而题目给出的角度、我们口头说的“30度”、“方位角120度”都是角度制。

  • 我们的教训:在代码中明确区分。所有从题目读取或由太阳位置算法计算出的角度,在代入三角函数计算前,必须用math.radians()转换为弧度。反之,如果需要输出角度,则用math.degrees()转换回来。我们在初期因为忘了转换,导致计算出的向量完全错误,效率值出现大于1或为负的荒谬结果。

5.3 截断效率模型的适用性与简化边界

我们采用的高斯近似解析公式虽然简洁,但它有几个强假设:光斑能量呈圆对称高斯分布、接收器为圆形、光学误差服从高斯分布。实际情况可能更复杂。

  • 应对策略
    1. 模型说明:在论文中必须明确指出这些假设,并讨论其合理性。例如,对于矩形接收器,该公式需要修正。
    2. 参数敏感性分析:在问题一求解后,可以简单分析截断效率对关键参数(如距离D、误差σ_slope)的敏感性。这能为后续优化问题提供洞察,例如,远离塔的镜子可能因截断损失过大而不经济。
    3. 备用方案:如果时间允许,可以提及更精确的蒙特卡洛光线追迹法是行业标准,但计算成本高,本模型采用解析近似以平衡精度与速度。

5.4 效率乘积模型的独立性假设

我们的总效率是四个子效率的乘积,这隐含了假设:这些效率因子是相互独立的。实际上,它们可能存在弱耦合。例如,阴影遮挡严重的位置,其余弦效率通常也较低(因为镜子可能处于边缘位置,指向角不佳)。但在问题一的瞬时计算中,这种独立性假设是普遍接受且合理的。

  • 在论文中的处理:明确指出这一假设,并说明其对于评估瞬时性能是可行的。在后续问题(如年化计算、布局优化)中,如果需要更精确的年度总能量评估,则需要考虑太阳位置变化下各效率因子的时间相关性。

6. 从问题一延伸:模型的价值与后续优化接口

完成问题一的建模与求解,绝不仅仅是算出一堆效率数字。它的真正价值在于为整个赛题搭建了一个坚实、可扩展的计算核心。

6.1 模型输出的深度利用

计算出的单镜效率矩阵,可以可视化出来,生成镜场效率分布云图。用matplotlibcontourfscatter(颜色映射效率值)可以直观显示:哪些区域的镜子效率高(可能靠近塔中心,距离适中,指向角好),哪些区域效率低(边缘、距离过远或过近)。这张图本身就是对镜场设计合理性的一个初步诊断。

6.2 为问题二、三奠定基础

问题一模型是一个“原子”模型。在问题二(计算特定时刻的镜场总输出功率)中,我们只需遍历所有镜子,将每面镜子的效率乘以镜面面积和该时刻的法向直接辐射辐照度(DNI),再求和即可。公式大致为:总功率 = DNI * 镜面总面积 * 镜场平均光学效率(更精确的是对各镜求和)。 在问题三(优化镜场布局以最大化年均输出)中,问题一的模型就成为了核心的目标函数计算器。优化算法(如遗传算法、粒子群算法)每提出一个新的镜场坐标布局方案,都需要调用问题一的模型来计算该布局在多个代表性时间点(如春分、夏至、秋分、冬至的多个小时)的效率,进而估算年总输出。此时,计算速度就至关重要,这也是为什么我们在问题一就采用解析模型而非光线追迹的原因。

6.3 模型的可扩展性思考

一个健壮的模型应该易于扩展。我们在编程时,将效率计算封装成了独立的函数,参数清晰。如果后续需要考虑阴影遮挡(η_shadow),我们可以编写一个calculate_shadow_efficiency(mirror_list, sun_vec)函数,判断每面镜子的受影情况,然后无缝集成到总效率计算中。同样,如果接收器是平面矩形而非圆形,截断效率函数也需要相应调整。这种模块化的设计思路,在三天紧张的竞赛中,能极大提高代码的可维护性和迭代效率。

回顾整个问题一的解决过程,从理解物理背景、拆解效率因子,到建立数学模型、编程实现,再到排查陷阱、思考延伸,这正是一个完整的数学建模实战闭环。它锻炼的不仅仅是数学和编程能力,更是将复杂工程问题抽象化、条理化的思维能力。希望这份详细的复盘,能为你理解这类优化设计问题提供一个扎实的起点。记住,好的开始是成功的一半,把问题一的模型做扎实了,后面的路会好走很多。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询