CN115049519A - 一种水上光伏对下覆开放水域环境水温影响的计算方法 - Google Patents

Abstract

本发明涉及一种水上光伏对下覆开放水域环境水温影响的计算方法，包括：收集信息；构建网格；设定计算时间参数；水动力边界条件输入；计算流速分布；水温分布计算准备；计算水温分布；判断是否达到采集边界。本发明通过水体热量平衡计算方法，以及一系列各通量分项影响因子受光伏阵列影响的计算公式，并以光伏阵列内外水域各项不同数据的比较，针对开放水域内的桩基水上光伏电站，依据当地的气象条件和即将建设的光伏阵列布设参数，定量计算水上光伏电站对水温变化的影响，根据此方法计算的水温及水温变化可作为评估光伏电站对水环境影响的依据，为将来光伏电站建成后，对该水域的开发利用提供了科学依据。

Description

一种水上光伏对下覆开放水域环境水温影响的计算方法

技术领域

本发明涉及一种水上光伏对下覆开放水域环境水温影响的计算方法，是一种水温观测计算方法，是一种用于水面光伏电站的光伏阵列影响的计算方法。

背景技术

水上光伏电站在水面上布设大量密集的光伏板阵列，光伏板阵列通过遮蔽和吸收太阳短波辐射、改变局部风场等方式，影响水面接受的辐射热和水表蒸发、对流等换热过程，进而改变水体温度。作为水生态的重要环境因子，水温是评估水上光伏电站对水环境影响的重要指标。目前主要基于有限的观测进行光伏电站对水温影响规律的定性分析，缺少系统性的定量计算和评价方法。

发明内容

为了克服现有技术的问题，本发明提出了一种水上光伏对下覆开放水域环境水温影响的计算方法。所述的方法针开放水域内的桩基水上光伏电站，提出了依据气象条件和光伏阵列布设参数，定量计算水上光伏电站中水温变化的方法，根据此方法计算的水温及水温变化可作为评估光伏电站对水环境影响的依据。

本发明的目的是这样实现的：一种水上光伏对下覆开放水域环境水温影响的计算方法，所述的方法的步骤如下：

步骤1，收集信息：收集拟建造的光伏阵列的各项设计数据和当地水文气象资料；光伏阵列的设计数据包括：光伏板尺寸、光伏板间距、支柱分布间距；当地水文气象资料包括：太阳辐射的总辐射量I，通过计算或收集获取的散射辐射量I_d；大气长波辐射、气温、云量、湿度、气压、水面风速、水温，以及地形资料；

当地水文气象资料是指在没有建造光伏阵列之前该区域的天气、水文资料和数据，或是在已经建成的光伏阵列之外的水域进行观测的天气、水文资料或数据，本实施例为区别，将没有光伏阵列或在光伏阵列之外收集和观测的资料或数据称之为“光伏区外水域”数值，代表该数值的字母不加脚标；将光伏阵列中收集和观测的资料或数据称之为“光伏区内水域”数值，代表该数值的字母加脚标pv，以方便识别。

步骤2，构建网格：对拟建造光伏阵列的受影响水域建立水动力数学模型网格；

建立数学模型网格是为了构建计算单元，在一个网格内进行温度计算，对于本实施例来说，一个网格可以设置为边长5米的正方形，投影在水面上，形成网格。

步骤3，设定计算时间参数：计算时间参数包括：计算起始时间t₀，时间步长，计算结束时间t_n；

步骤4，水动力边界条件输入：水动力边界条件包括：水边界、底部摩擦项；

子步骤4.1，确定水边界：采用输入潮位边界，采用TPXO全球海潮模型给出水边界上的潮汐调和常数作为输入；

子步骤4.2，确定底部摩擦项：

i.光伏区外水域：设定Chezy系数C_2D作为底部摩擦输入项，床面剪切力τ_b：

式中：U为水深平均水平向流速；ρ₀为水的密度；

ii.光伏区内水域：给定Chezy系数C_2D以及单位面积上的桩密度N_a，桩柱高度h，桩直径d，以及桩群的平均阻力系数C_D值；

对光伏区水域，床底总的切应力τ′_b：

步骤5，计算流速分布：从t₀时刻开始，求解水动力方程；所述的水动力方程包括：水深平均的连续方程和模型动量方程；

水深平均的连续方程：

其中：t为当前时间；ξ和η为曲线坐标系下水平方向的坐标；G_ξξ、G_ηη为曲线坐标系转换为直角坐标系的转换系数；U、V分别为ξ和η方向的垂向平均流速；Q表示源汇项；

模型动量方程：

式中：u、v、ω是曲线正交坐标系下的流速分量；f是Coriolis参数；F_ξ和F_η分别为ξ和η方向的紊动动量通量；P_ξ和P_η为ξ和η方向的水压力梯度；ν_v为垂向粘滞系数；M_ξ和M_η为ξ和η方向上动量的源汇项；

水动力模型基于Boussinesq和垂向静水压力假定，采用三维曲线坐标系，垂向坐标表示如下：

其中：d为参考平面以下的水深值，H是全水深，z为物理空间垂向坐标，ζ为自由表面高程；σ_h为相对水深；

步骤6，水温分布计算准备：水温分布计算准备包括：水边界给定已知水温；底边边界给定热通量为0；给定自由表面热通量边界条件；

所述的自由表面热通量边界条件包括：太阳短波辐射、大气长辐射、蒸发散射、对流热通量，气温，湿度，风速，大气压力；通过光伏遮蔽效应求解得到光伏区的自由表面热通量：

6.1太阳短波辐射的计算：

i.光伏区外水域的太阳短波辐射采用总辐射量I的实测值；

ii.光伏区内水域的太阳短波辐射I_pv计算包括如下子步骤：

子步骤6.1.1，根据收集的信息整理获得太阳辐射总辐射量I，并计算散射辐射量I_d；

子步骤6.1.2，计算天空视域比例因子CF_sky：

式中θ₁、θ₂、θ₃、θ₄分别为计算空隙后板和前板下檐、上檐到积分点的夹角；n为积分区间L的分段数，是阵列前后两排光伏板之间从前排前檐到后排前檐；ml_max为积分区间左侧最大空隙数；mr_max为积分区间右侧含区分区间本身的最大空隙数；

子步骤6.1.3，计算光伏板后檐对应阴影宽度D_h：

D_h＝H_pv·cotα·cosγ

式中：α为太阳高度角，γ为太阳方位角，H_pv为光伏板上下沿高度差；

子步骤6.1.4，计算透过比例因子CF_dir：直射辐射在前后阵列之间宽度为L内的透过比例因子CF_dir为：

其中：D_pv为光伏板前檐到后檐在水平面投影面中的距离；

子步骤6.1.5，计算光伏区下垫面总辐射量I_pv：

I_pv＝(I-I_d)CF_dir+I_dCF_sky；

6.2大气长辐射的计算：

i.光伏区外水域的大气长波辐射采用实测值或根据实测气温、云量计算；

光伏区外的水域的大气长波辐射E_a计算公式：

E_a＝εσT_a0 ⁴g(Fc)

式中：ε为发射率；T_a0为t₀时刻气温；g(Fc)为云量Fc的函数：

g(Fc)＝(1+0.17Fc²)；

ii.光伏区内水域的大气长波辐射E_apv，根据实测或计算得到；

计算光伏区内水域的大气长波辐射包括如下子步骤：

子步骤6.2.1，光伏区内水域等效云量系数计算：根据光伏阵列对天空遮挡的视域因子，以及收集的云量信息求解光伏区内水域考虑光伏板遮挡后的等效云量系数Fc′：

Fc′＝(1-CF_sky)+CF_sky·Fc；

子步骤6.2.2：计算有光伏影响下的等效云量函数：

g(Fc′)＝(1+0.17Fc′²)；

子步骤6.2.3：计算有光伏影响下的大气长波辐射E_apv：

E_apv＝E_a·g(F′_C)/g(F_c)；

6.3蒸发散射的计算：

i.光伏区外水域的蒸发散热计算：

子步骤6.3.1，计算水面饱和水汽压力、空气中水汽压力：

水面饱和水汽压力e_s：

e_s＝f(T_s)

其中：e_a空气中水汽压力：

e_a＝f(T_a)

T_s为水表面温度；T_a为气温；

水汽压力f(T)计算公式：

其中：T为水表面温度或气温；E₀为纯水平液面饱和水汽压，计算公式：

其中：T_c为水的三相点温度；T为水表面温度或气温；

子步骤6.3.2，计算蒸发量：

i.对于光伏区外水域，蒸发公式：

E＝f(u)·(e_s-e_a)

式中：f(u)为风速函数，u为光伏区外水域的水面风速，风速函数f(u)：

ii.对于光伏区内水域的水面风速u_pv计算公式：

u_pv＝r_s·u

其中：r_s＝u₀/u_c0；

u₀：已建成的类似光伏阵列内地面或水面风速；u_c0：已建成的类似光伏阵列外地面或水面风速；

光伏区内水域蒸发量：

E_pv＝f(u_pv)·(e_s-e_a)；

子步骤6.3.3，计算蒸发潜热：蒸发潜热L_v计算公式：

L_v＝2.5×10⁶-2.3×10³·T_s0；

子步骤6.3.4，蒸发潜热通量Q_E通过单位时间蒸发量得到：

i.光伏区外水域：Q_E＝L_vρ₀E

ii.光伏区内水域：Q_Epv＝L_vρ₀E_pv

其中：ρ₀为水体密度(kg/m³)

6.4对流热通量的计算：

子步骤6.4.1：计算波文比b：

b＝0.62P_a/1000(T_s－T_a)/(e_s－e_a)

其中：P_a为大气压力；

子步骤6.4.2：根据波文比计算对流热通量：

i.光伏区外水域：Q_H＝b·Q_E；

ii.光伏区内水域：Q_Hpv＝b·Q_Epv；

6.5水体向外长波辐射热通量E_W：水体向外长波辐射热通量计算公式：

E_w＝εσT_s ⁴；

6.6计算t时刻自由表面总热通量F：

i.光伏区外水域：

F＝I(1-α_s)+E_a-E_w-Q_E-Q_H；

ii.光伏区内水域：

F_pv＝I_pv(1-α_s)+W_apv-W_w-Q_Epv-Q_Hpv

其中，α_s为太阳短波辐射反射比；

步骤7，计算水温分布：求解物质运输方程，获得t时刻水温分布：

求解物质输运方程：

式中：c为物质浓度；D_H为水平扩散系数；D_v为垂向扩散系数；λ_d为物质衰变系数；S为源项；

步骤8，判断是否达到采集边界：判断是否t＝t_n，如果是，则结束计算，如果否则回到步骤5计算下一时段的水温分布。

进一步的，子步骤6.1.1所述的环境散射辐射量的计算方式包括如下子步骤：

子步骤6.1.1.1，计算大气层外辐射量I₀：

其中：E_sc为太阳常数；δ为赤纬角；

为当地纬度；ω为时角；λ为日-地距离变化引起大气层上界太阳辐照度的修正值：

其中：n_d为一年中的日序号；

子步骤6.1.1.2，计算晴空指数：

子步骤6.1.1.3，计算散射辐射与总辐射之比：下垫面散射辐射I_d与总辐射I之比的计算公式：

子步骤6.1.1.4，计算散射辐射：

I_d＝I·r_d。

本发明的优点和有益效果是：本发明通过水体热量平衡计算方法，以及一系列各通量分项影响因子受光伏阵列影响的计算公式，并以光伏阵列内外水域各项不同数据的比较，针对开放水域内的桩基水上光伏电站，依据当地的气象条件和即将建设的光伏阵列布设参数，定量计算水上光伏电站对水温变化的影响，根据此方法计算的水温及水温变化可作为评估光伏电站对水环境影响的依据，为将来光伏电站建成后，对该水域的开发利用提供了科学依据。

附图说明

下面结合附图和实施例对本发明作进一步说明。

图1是本发明实施例一所述方法多块光伏板公式的参数标注示意图；

图2是本发明实施例一所述方法的流程图；

图3是本发明实施例二所述应用实例的光伏阵列的平面布置示意图；

图4是本发明实施例二所述应用实例的光伏对底部阻力X方向糙率值影响结果示意图；

图5是本发明实施例二所述应用实例的光伏对底部阻力Y方向糙率值影响结果示意图；

图6是本发明实施例二所述应用实例计算统计的增加光伏后自然海域水温的降温等值线分布。

具体实施方式

实施例一：

本实施例是一种水上光伏对下覆开放水域环境水温影响的计算方法。

开放水域存在水流运动，流动会使得热量在不同空间内进行交换，因此建造于开放水域的光伏电站，其对水温的影响不局限在光伏所在的下覆水域，也会影响到光伏区外围的水域，这是影响范围与封闭水域的不同。开放水域光伏区水体受光伏遮蔽影响，表面热通量项中的太阳辐射、大气辐射会发生变化，光伏对风的遮挡会导致风速发生变化，进而表面热通量项中的蒸发和对流热通量发生变化，这是光伏造成的水温边界条件的变化，本实施例对光伏的水温边界条件通过光伏板对散射辐射和直射辐射的折减比例求解纳入光伏造成的太阳总辐射变化、将光伏板视为全遮挡纳入等效云量系数反映光伏区大气长波辐射的变化，通过风速折减系数求解光伏区蒸发和对流散热量的变化。光伏板密集的支撑柱会对水流运动产生影响，进而影响到热量的输运，单纯支撑柱本身会对开放水域原有的水温分布结构带来变化，支撑柱对水流运动的影响通过光伏造成的水动力底部摩擦边界条件变化输入到水动力求解过程中。本实施例将桩柱阻力叠加到底部摩擦阻力中，以反映桩柱对水流运动的影响。计算要点如下：

一、开放水域水动力和水温模拟计算方法：

水域的水动力求解一般基于N-S方程进行数值求解，求解的控制方程包括连续方程和动量方程、紊动方程，水温采用物质输运方程求解。以Delft3D软件为例，水动力模型基于Boussinesq和垂向静水压力假定，采用三维曲线坐标系，垂向坐标表示如下：

其中d为参考平面以下的水深值，H是全水深，z为物理空间垂向坐标，ζ为自由表面高程(参考平面以上的数值)。在底层σ_h＝-1，表层σ_h＝0。

水深平均的连续方程：

其中，t为时间；ξ和η为曲线坐标系下水平方向的坐标；G_ξξ、G_ηη为曲线坐标系转换为直角坐标系的转换系数；U、V分别为ξ和η方向的垂向平均流速；Q表示源汇项。

模型动量方程：

式中u、v、ω是曲线正交坐标系下的流速分量；f是Coriolis参数；F_ξ和F_η分别为ξ和η方向的紊动动量通量；P_ξ和P_η为ξ和η方向的水压力梯度，ν_v为垂向粘滞系数，M_ξ和M_η为ξ和η方向上动量的源汇项。

模拟中采用斜压模型，水体密度为温度和盐度的函数，采用UNESCO状态方程表示，每个求解步通过状态方程求解温、盐度分布引起的密度差，在压力项P中予以反映。

物质输运方程：

式中c为物质浓度；D_H为水平扩散系数；D_v为垂向扩散系数；物质扩散系数的确定方法类比水流粘滞系数；λ_d为物质衰变系数(/day)；S为源项。

垂向紊动粘滞系数可通过紊流模型求解，Delft3D提供了Prandtl混合长模型，k-L模型和k-ε模型，本次计算选取k-ε模型。

在求解上述方程时需要给定边界条件，其中受光伏影响的包括求解水动力方程所需的底部阻力系数，求解物质输运方程需要给定水体表面热通量。

以Delft3D中提供的表面热交换模型(Murakami模型)为例，计算表面热通量各分项计算公式如下：

1)辐射热：

包括太阳辐射、大气辐射、水体向外辐射。

水体吸收的太阳辐射热φ_sn采用已知输入值。

大气辐射和水体向外辐射均为长波辐射，计算中将二者合为长波辐射φeb考虑，计算公式为：

其中：φ_w为水体对外辐射散热量；T_s为表面水温；T_a为气温；ε为水体透明系数，取值0.9；σ为Stefan-Boltzmann常数，σ＝5.67×10^-8(J/m²sK⁴)；e_a为空气中的水汽压力(hPa)；F_c为云量系数。

2)蒸发与对流散热：

蒸发散热φ_ev计算公式为：

φ_ev＝L_vρ₀E (7)

其中：ρ₀为水体密度，L_v为蒸发潜热，其单位为J/kg,表达式为：

L_v＝2.5·10⁶-2.3·10³T_s (8)

E为体积蒸发量：

E＝f(U₁₀)(e_s-e_a) (9)

e_s为与T_s对应的饱和水汽压力(hPa)，f(U₁₀)风速函数：

f(U₁₀)＝c_eU₁₀ (10)

U₁₀为水面上10m高处风速，m/s；c_e为Dalton系数，取值0.0012。

与之对应的对流散热计算公式为：

φ_c＝bφ_ev (11)

其中φ_c为对流散热量，b为波温比：

b＝0.66(T_s-T_a)/(e_s-e_a) (12)

3)自由表面热量交换总量：

φ_n＝φ_sn-φ_eb-φ_ev-φ_c (13)

式(13)中所述各热通量项受光伏遮挡均会发生变化，其中辐射热遮挡与板的空间布置相关，太阳辐射遮挡还会与太阳角度有关。蒸发和对流散热计算包含了气温、湿度和风速，其中受光伏板影响较为显著的是风速。因此受光伏影响的水温变化主要考虑太阳辐射、大气长波辐射的遮挡效应以及风速的遮挡效应进行计算。

二、表面热通量受光伏影响的计算方法：

各分项影响因子受光伏影响的计算公式：

1.太阳辐射：

散射辐射为各向同性漫射，折减系数采用天空视域因子，对多排缝隙判断该缝隙的空域是否入射到计算单元以及计算入射角度的方法如下：

特定缝隙的散射光能否进入积分单元要判断缝隙左侧板的最低入射角是否大于右侧板的最高入射角，如果是,则该二夹角间的空域为积分单元上可视空域，光线可达，否则该区域空隙被前后板互相遮挡，无法到达积分单元，总的天空视域因子计算公式如下：

式中：θ₁、θ₂、θ₃、θ₄分别为计算空隙后板和前板下檐、上檐到积分点的夹角，如图1所示；n为积分区间L的分段数；ml_max为积分区间左侧最大空隙数；mr_max为积分区间右侧(含区分区间本身)最大空隙数。图1显示了光伏阵列某个特定单元的垂向立面图，简略的显示了四排光伏板，并将中间的两个光伏板按照由南向北的方向设定为前排光伏板1和后排光伏板2，并将阳光入射方向3设定在入射角α的角度上，将积分点4设在前排光伏板上檐在水面5的投影点(线)上。

直射辐射入射比例可根据太阳高度角和太阳方位角求解光照区比例获得。对于特定单元，当单元位于阴影区时，直射辐射透过比为0，当单元位于光照区时，直射辐射透过比为1。对于整个光伏区，前后两阵列之间光照区宽度占间距的比值即为直射透过比例。

光伏板后檐对应阴影宽度：

D_h＝H_pv·cotα·cosβ (14)

式中α为太阳高度角，β为太阳方位角，与时间有关，可采用相关公式求解。

直射辐射在前后阵列宽度L内的透过比例为：

光伏区根据输入的总辐射I，利用地面散射辐射与总辐射之比的计算公式得到散射辐射值I_d，

根据输入的总辐射I和散射辐射值I_d，按照式(13)和(15)计算方法分别计算光伏区散射辐射和直射辐射折减系数，折减后的值相加得到光伏区总辐射。

I_pv＝(I-I_d)CF_dir+I_dCF_sky (17)

2.长波辐射：

光伏遮挡会减少大气向水体的辐射，大气辐射公式的一般形式为：

E_a＝εσT_a ⁴g(Fc) (18)

式中ε为发射率，可通过实测气象资料进行率定；σ为玻尔兹曼常数，取值5.67×10^-8；T_a为气温，g(Fc)为云量Fc的函数，常用的形式为：

g(Fc)＝(1+0.17Fc²) (19)

光伏对天空遮挡的比例系数F_pv可以采用太阳短波散射辐射计算的天空视域因子CF_sky得：

F_pv＝1-CF_sky (20)

将光伏遮挡区视为云量，光伏区考虑光伏板遮挡后的云量系数为：

Fc′＝F_pv+CF_sky·F_c (21)

将Fc′带入(19)式，再按(18)式求解得到光伏区受遮挡的大气辐射。

3风速：

光伏板对风的遮挡与光伏板布置间距、倾角、高度，以及风向等均相关，可以通过原型观测统计分析，或者模型试验或数值计算方法获取折减系数。之后带入(10)式求解光伏遮挡区的蒸发和对流散热。

三、光伏支撑柱对水流运动的阻力计算：

对全水深平均的二维流动，床面剪切力采用二次摩擦定律：

式中：U为水深平均水平向流速；ρ₀为水的密度；C_2D为二维Chezy系数(m^1/2/s)。

桩群阻力：

式中：C_D为桩群的平均阻力系数。

网格尺寸为a×a范围内的桩群，沿水流方向的阻力切应力为：

τ_z＝F/a² (24)

床底总的切应力在原有基础上增加了桩柱阻力切应力，变为：

在动量方程中将上式带入计算阻力，即可对光伏区的阻力系数进行修正。

对三维模型，床面切应力与近床底层流速相关：

式中：u_b为近底层水平流速；C_3D为阻力系数。可采用与二维同样的方法对阻力系数进行修正。

本实施例所述方法的具体步骤如下：

步骤1，收集信息：收集拟建造的光伏阵列的各项设计数据和当地水文气象资料；光伏阵列的设计数据包括：光伏板尺寸、光伏板间距、支柱分布间距；当地水文气象资料包括：太阳辐射的总辐射量I，通过计算或收集获取的散射辐射量I_d；大气长波辐射、气温、云量、湿度、气压、水面风速、水温，以及地形资料。

当地水文气象资料是指在没有建造光伏阵列之前该区域的天气、水文资料和数据，或是在已经建成的光伏阵列之外的水域进行观测的天气、水文资料或数据，本实施例为区别，将没有光伏阵列或在光伏阵列之外收集和观测的资料或数据称之为“光伏区外水域”数值，代表该数值的字母不加脚标；将光伏阵列中收集和观测的资料或数据称之为“光伏区内水域”数值，代表该数值的字母加脚标pv，以方便识别。

步骤2，构建网格：对拟建造光伏阵列的受影响水域建立水动力数学模型网格。

建立数学模型网格是为了构建计算单元，在一个网格内进行温度计算。开放水域，例如近海海域，水温分布的实时变化是大范围气象和水文场耦合作用的结果，受制于外海水温资料精确获取的困难，为了尽可能减小开边界水温条件赋值对工程海域温度模拟结果的影响，所构建的网格模拟边界应尽量远离工程区。例如：可以将水温模拟范围在拟建造光伏电站选址处约250km范围内，模拟区域平面采用4重嵌套网格，最内层网格尺度50m，最外层网格尺度约2000m，平面网格布置见图3所示，垂向上采用6层σ网格。

步骤3，设定计算时间参数：计算时间参数包括：计算起始时间t₀，时间步长，计算结束时间t_n。

通常可以将计算起始时间设定为每天早晨太阳升起前，例如早上5点钟，可以设定时间步长为半小时，计算结束时间为每天下午太阳落山时，例如下午6点。

步骤4，水动力边界条件输入：水动力边界条件包括：水边界、底部摩擦项。

子步骤4.1，确定水边界：采用输入潮位边界，采用TPXO全球海潮模型给出水边界上的潮汐调和常数作为输入。

子步骤4.2，确定底部摩擦项：

i.光伏区外水域：设定Chezy系数(m^1/2/s)C_2D作为底部摩擦输入项，床面剪切力τ_b：

式中：U为水深平均水平向流速(计算过程中计算得到的值)；ρ₀为水的密度。

ii.光伏区内水域：给定Chezy系数(m^1/2/s)C_2D以及单位面积上的桩密度N_a，桩柱高度h，桩直径d，以及桩群的平均阻力系数C_D值。

对光伏区水域，床底总的切应力τ′_b在原有床面切应力基础上增加桩柱阻力切应力，变为：

步骤5，计算流速分布：从t₀时刻开始，求解水动力方程；所述的水动力方程包括：水深平均的连续方程和模型动量方程。

水深平均的连续方程：

其中：t为当前时间；ξ和η为曲线坐标系下水平方向的坐标；G_ξξ、G_ηη为曲线坐标系转换为直角坐标系的转换系数；U、V分别为ξ和η方向的垂向平均流速；Q表示源汇项。

模型动量方程：

式中：u、v、ω是曲线正交坐标系下的流速分量；f是Coriolis参数；F_ξ和F_η分别为ξ和η方向的紊动动量通量；P_ξ和P_η为ξ和η方向的水压力梯度；ν_v为垂向粘滞系数；M_ξ和M_η为ξ和η方向上动量的源汇项。

水域的水动力求解一般基于N-S方程进行数值求解，求解的控制方程包括连续方程和动量方程、紊动方程，水温采用物质输运方程求解。以Delft3D软件为例，水动力模型基于Boussinesq和垂向静水压力假定，采用三维曲线坐标系，垂向坐标表示如下：

其中：d为参考平面以下的水深值，H是全水深，z为物理空间垂向坐标，ζ为自由表面高程(参考平面以上的数值)；σ_h为相对水深；在底层σ_h＝-1，表层σ_h＝0。

步骤6，水温分布计算准备：水温分布计算准备包括：水边界给定已知水温；底边边界给定热通量为0；给定自由表面热通量边界条件；

所述的自由表面热通量边界条件包括：太阳短波辐射、大气长辐射、蒸发散射、对流热通量，气温，湿度，风速，大气压力；通过光伏遮蔽效应求解得到光伏区的自由表面热通量：

6.1太阳短波辐射的计算：

i.光伏区外水域的太阳短波辐射采用总辐射量I的实测值；

ii.光伏区内水域的太阳短波辐射I_pv采用无光伏影响的实测值和计算的折减系数得到。

光伏区内水域的太阳短波辐射I_pv计算包括如下子步骤：

子步骤6.1.1，根据收集的信息整理获得太阳辐射总辐射量I，并计算散射辐射量I_d；

子步骤6.1.2，计算天空视域比例因子CF_sky：

式中θ₁、θ₂、θ₃、θ₄分别为计算空隙后板和前板下檐、上檐到积分点的夹角；n为积分区间L的分段数；ml_max为积分区间左侧最大空隙数；mr_max为积分区间右侧含区分区间本身的最大空隙数。

子步骤6.1.3，计算光伏板后檐对应阴影宽度D_h：

D_h＝H_pv·cotα·cosγ

式中：α为太阳高度角，γ为太阳方位角，H_pv为光伏板上下沿高度差。

子步骤6.1.4，计算透过比例因子CF_dir：直射辐射在前后阵列之间宽度为L内的透过比例因子CF_dir为：

其中：L为：阵列前后两排光伏板之间从前排前檐到后排前檐；D_pv为光伏板前檐到后檐在水平面投影面中的距离。

子步骤6.1.5，计算光伏区下垫面总辐射量I_pv：

I_pv＝(I-I_d)CF_dir+I_dCF_sky。

6.2大气长辐射的计算：

i.光伏区外水域的大气长波辐射采用实测值或根据实测气温、云量计算；

根据收集的信息获得实测大气长波辐射或通过实测的气温、云量计算获得无光伏影响的大气长波辐射值。光伏区外的水域的大气长波辐射E_a计算公式：

E_a＝εσT_a0 ⁴g(Fc)

式中：ε为发射率，可通过实测气象资料进行率定，默认值取0.97，σ为玻尔兹曼常数，取值5.67×10^-8(W˙m^-2˙℃^-4)；T_a0为t₀时刻气温(K)；g(Fc)为云量Fc的函数：

g(Fc)＝(1+0.17Fc²)。

ii.光伏区内水域的大气长波辐射E_apv，根据实测或计算得到，光伏区内水域的大气长波辐射E_apv比光伏区外水域的大气长波辐射E_a多了求解等效云量系数。

计算光伏区内水域的大气长波辐射包括如下子步骤：

子步骤6.2.1，光伏区内水域等效云量系数计算：根据光伏阵列对天空遮挡的视域因子，以及收集的云量信息求解光伏区内水域考虑光伏板遮挡后的等效云量系数Fc′：

Fc′＝(1-CF_sky)+CF_sky·Fc；

子步骤6.2.2：计算有光伏影响下的等效云量函数：

g(Fc′)＝(1+0.17Fc′²)；

子步骤6.2.3：计算有光伏影响下的大气长波辐射E_apv：

E_apv＝E_a·g(F′_C)/g(F_c)；

6.3蒸发散射的计算：

i.光伏区外水域的蒸发散热计算：

子步骤6.3.1，计算水面饱和水汽压力、空气中水汽压力：

水面饱和水汽压力e_s：

e_s＝f(T_s)

其中：e_a空气中水汽压力：

e_a＝f(T_a)

T_s为水表面温度(K)；T_a为气温(K)；

水汽压力f(T)计算公式：

其中：T为水表面温度或气温(K)；E₀为纯水平液面饱和水汽压(hPa)，计算公式：

其中：T_c为水的三相点温度，T_c＝273.16K；T为水表面温度或气温(K)；

子步骤6.3.2，计算蒸发量：

i.对于光伏区外水域，选用经典的道尔顿蒸发公式形式：

E＝f(u)·(e_s-e_a)

式中：f(u)为风速函数，u为光伏区外水域的水面风速(m/s)，即没有受到光伏阵列影响的风速，风速函数f(u)可以采用国内常用的闵骞修正的风速函数：

ii.对于光伏区内水域的水面风速u_pv在无光伏影响的风速u基础上计入光伏遮挡影响，u_pv据类似已建光伏区实测风速折减比r_s和拟计算的无光伏影响风速u计算，即找一处或几处(包括地面或水面)已经建成的光伏阵列，检测光伏阵列内、外的风速，用以计算风速递减比r_s，计算公式：

u_pv＝r_s·u

其中：r_s＝u₀/u_c0；

u₀：已建成的类似光伏阵列内地面或水面风速；u_c0：已建成的类似光伏阵列外地面或水面风速；

光伏区内水域蒸发量：

E_pv＝f(u_pv)·(e_s-e_a)；

子步骤6.3.3，计算蒸发潜热：蒸发潜热L_v(J/kg)计算公式：

L_v＝2.5×10⁶-2.3×10³·T_s0；

子步骤6.3.4，蒸发潜热通量Q_E通过单位时间蒸发量得到：

i.光伏区外水域：Q_E＝L_vρ₀E

ii.光伏区内水域：Q_Epv＝L_vρ₀E_pv

其中：ρ₀为水体密度(kg/m³)

6.4对流热通量的计算：

子步骤6.4.1：计算波文比b：

b＝0.62P_a/1000(T_s－T_a)/(e_s－e_a)

其中：P_a为大气压力(hPa)；

子步骤6.4.2：根据波文比计算对流热通量：

i.光伏区外水域：Q_H＝b·Q_E；

ii.光伏区内水域：Q_Hpv＝b·Q_Epv；

6.5水体向外长波辐射热通量E_W：水体向外长波辐射热通量计算公式：

E_w＝εσT_s ⁴；

6.6计算t时刻自由表面总热通量F：

i.光伏区外水域：

F＝I(1-α_s)+E_a-E_w-Q_E-Q_H；

ii.光伏区内水域：

F_pv＝I_pv(1-α_s)+E_apv-E_w-Q_Epv-Q_Hpv

其中，α_s为太阳短波辐射反射比，可通过实测资料率定，默认取值0.06。

步骤7，计算水温分布：求解物质运输方程，获得t时刻水温分布：

求解物质输运方程：

式中：c为物质浓度；D_H为水平扩散系数；D_v为垂向扩散系数；物质扩散系数的确定方法类比水流粘滞系数。λ_d为物质衰变系数(/day)；S为源项。

垂向紊动粘滞系数可通过紊流模型求解，Delft3D提供了Prandtl混合长模型，k-L模型和k-ε模型，本实施例计算选取k-ε模型。

步骤8，判断是否达到采集边界：判断是否t＝t_n，如果是，则结束计算，如果否则回到步骤5计算下一时段的水温分布。

本步骤判断是否达到采集边界，例如：以天为计算边界的，如果设定下午6点，太阳已经落山，太阳辐射为零时为采集边界，达到这一边界则结束计算。

实施例二：

本实施例是实施例一的改进，是实施例一关于环境辐射值计算方法的细化。所述子步骤6.1.1所述的环境散射辐射量的计算方式包括如下子步骤：

子步骤6.1.1.1，计算大气层外辐射量I₀：

其中：E_sc为太阳常数；δ为赤纬角；

为当地纬度；ω为时角；λ为日-地距离变化引起大气层上界太阳辐照度的修正值：

其中：n_d为一年中的日序号；

子步骤6.1.1.2，计算晴空指数：

子步骤6.1.1.3，计算散射辐射与总辐射之比：下垫面散射辐射I_d与总辐射I之比的计算公式：

子步骤6.1.1.4，计算散射辐射：

I_d＝I·r_d。

应用实例：

实施案例位于江苏省连云港市拟建项目滩涂海域，光伏阵列的布置形式如图3所示。光伏阵列的组串支架布置形式按竖向2行28列(图3中仅画出2列)，均采用最佳倾角28°固定安装在支架上。支架桩柱8采用12m长直径400mm的PHC-AB管桩，桩柱间南北向间距9.0m，东西向间距为3.8m，两个相邻光伏板的第一个桩柱间距为2.364m，光伏板最低点标高5.70m。

选取的典型月实际逐小时过程，潮汐条件采用实际潮汐过程，计算一个月，前半月作为预热，选取后半个月半月潮作为结果输出时段。计算的光伏桩柱对水流底部阻力的影响结果见图4、5所示。其中图4是光伏对底部阻力X方向糙率值影响结果示意图，图是光伏对底部阻力Y方向糙率值影响结果示意图。有图4、5可以看出南北向阻力相比东下向增加较为明显，这主要是由于光伏阵列朝南布置，阵列内东西向桩柱间距仅4m，以致南北向迎流面桩柱密度更大。

图6为计算统计的增加光伏后自然海域水温的降温等值线分布，靠近湾顶平均降温幅度最大，可能1.2℃，向外海侧降温幅度逐渐下降，0.2℃温降影响范围大于光伏区布置范围。

最后应说明的是，以上仅用以说明本发明的技术方案而非限制，尽管参照较佳布置方案对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案(比如光伏阵列的排布方式、各种公式的运用、步骤的先后顺序等)进行修改或者等同替换，而不脱离本发明技术方案的精神和范围。

Claims (2)

1.一种水上光伏对下覆开放水域环境水温影响的计算方法，其特征在于，所述的方法的步骤如下：

步骤1，收集信息：收集拟建造的光伏阵列的各项设计数据和当地水文气象资料；光伏阵列的设计数据包括：光伏板尺寸、光伏板间距、支柱分布间距；当地水文气象资料包括：太阳辐射的总辐射量I，通过计算或收集获取的散射辐射量I_d；大气长波辐射、气温、云量、湿度、气压、水面风速、水温，以及地形资料；

当地水文气象资料是指在没有建造光伏阵列之前该区域的天气、水文资料和数据，或是在已经建成的光伏阵列之外的水域进行观测的天气、水文资料或数据，本实施例为区别，将没有光伏阵列或在光伏阵列之外收集和观测的资料或数据称之为“光伏区外水域”数值，代表该数值的字母不加脚标；将光伏阵列中收集和观测的资料或数据称之为“光伏区内水域”数值，代表该数值的字母加脚标pv，以方便识别。

步骤2，构建网格：对拟建造光伏阵列的受影响水域建立水动力数学模型网格；

建立数学模型网格是为了构建计算单元，在一个网格内进行温度计算，对于本实施例来说，一个网格可以设置为边长5米的正方形，投影在水面上，形成网格。

步骤3，设定计算时间参数：计算时间参数包括：计算起始时间T₀，时间步长，计算结束时间T_n；

步骤4，水动力边界条件输入：水动力边界条件包括：水边界、底部摩擦项；

子步骤4.1，确定水边界：采用输入潮位边界，采用TPXO全球海潮模型给出水边界上的潮汐调和常数作为输入；

子步骤4.2，确定底部摩擦项：

i.光伏区外水域：设定Chezy系数C_2D作为底部摩擦输入项，床面剪切力τ_b：

式中：U为水深平均水平向流速；ρ₀为水的密度；

ii.光伏区内水域：给定Chezy系数C_2D以及单位面积上的桩密度N_a，桩柱高度h，桩直径d，以及桩群的平均阻力系数C_D值；

对光伏区水域，床底总的切应力

步骤5，计算流速分布：从T₀时刻开始，求解水动力方程；所述的水动力方程包括：水深平均的连续方程和模型动量方程；

水深平均的连续方程：

其中：t为当前时间；ξ和η为曲线坐标系下水平方向的坐标；G_ξξ、G_ηη为曲线坐标系转换为直角坐标系的转换系数；U、V分别为ξ和η方向的垂向平均流速；Q表示源汇项；

模型动量方程：

式中：u、v、ω是曲线正交坐标系下的流速分量；f是Coriolis参数；F_ξ和F_η分别为ξ和η方向的紊动动量通量；P_ξ和P_η为ξ和η方向的水压力梯度；ν_v为垂向粘滞系数；M_ξ和M_η为ξ和η方向上动量的源汇项；

水动力模型基于Boussinesq和垂向静水压力假定，采用三维曲线坐标系，垂向坐标表示如下：

其中：d为参考平面以下的水深值，H是全水深，z为物理空间垂向坐标，ζ为自由表面高程；σ_h为相对水深；

步骤6，水温分布计算准备：水温分布计算准备包括：水边界给定已知水温；底边边界给定热通量为0；给定自由表面热通量边界条件；

所述的自由表面热通量边界条件包括：太阳短波辐射、大气长辐射、蒸发散射、对流热通量，气温，湿度，风速，大气压力；通过光伏遮蔽效应求解得到光伏区的自由表面热通量：

6.1太阳短波辐射的计算：

i.光伏区外水域的太阳短波辐射采用总辐射量I的实测值；

ii.光伏区内水域的太阳短波辐射I_pv计算包括如下子步骤：

子步骤6.1.1，根据收集的信息整理获得太阳辐射总辐射量I，并计算散射辐射量I_d；

子步骤6.1.2，计算天空视域比例因子CF_sky：

式中θ₁、θ₂、θ₃、θ₄分别为计算空隙后板和前板下檐、上檐到积分点的夹角；n为积分区间L的分段数，是阵列前后两排光伏板之间从前排前檐到后排前檐；ml_max为积分区间左侧最大空隙数；mr_max为积分区间右侧含区分区间本身的最大空隙数；

子步骤6.1.3，计算光伏板后檐对应阴影宽度D_h：

D_h＝H_pv·cotα·cosγ

式中：α为太阳高度角，γ为太阳方位角，H_pv为光伏板上下沿高度差；

子步骤6.1.4，计算透过比例因子CF_dir：直射辐射在前后阵列之间宽度为L内的透过比例因子CF_dir为：

其中：D_pv为光伏板前檐到后檐在水平面投影面中的距离；

子步骤6.1.5，计算光伏区下垫面总辐射量I_pv：

I_pv＝(I-I_d)CF_dir+I_dCF_sky；

6.2大气长辐射的计算：

i.光伏区外水域的大气长波辐射采用实测值或根据实测气温、云量计算；

光伏区外的水域的大气长波辐射E_a计算公式：

E_a＝εσT_a0 ⁴g(Fc)

式中：ε为发射率；T_a0为t₀时刻气温；g(Fc)为云量Fc的函数：

g(Fc)＝(1+0.17Fc²)；

ii.光伏区内水域的大气长波辐射E_apv，根据实测或计算得到；

计算光伏区内水域的大气长波辐射包括如下子步骤：

子步骤6.2.1，光伏区内水域等效云量系数计算：根据光伏阵列对天空遮挡的视域因子，以及收集的云量信息求解光伏区内水域考虑光伏板遮挡后的等效云量系数Fc′：

Fc′＝(1-CF_sky)+CF_sky·Fc；

子步骤6.2.2：计算有光伏影响下的等效云量函数：

g(Fc′)＝(1+0.17Fc′²)；

子步骤6.2.3：计算有光伏影响下的大气长波辐射E_apv：

E_apv＝E_a·g(F′_C)/g(F_c)；

6.3蒸发散射的计算：

i.光伏区外水域的蒸发散热计算：

子步骤6.3.1，计算水面饱和水汽压力、空气中水汽压力：

水面饱和水汽压力e_s：

e_s＝f(T_s)

其中：e_a空气中水汽压力：

e_a＝f(T_a)

T_s为水表面温度；T_a为气温；

水汽压力f(T)计算公式：

其中：T为水表面温度或气温；E₀为纯水平液面饱和水汽压，计算公式：

其中：T_c为水的三相点温度；T为水表面温度或气温；

子步骤6.3.2，计算蒸发量：

i.对于光伏区外水域，蒸发公式：

E＝f(u)·(e_s-e_a)

式中：f(u)为风速函数，u为光伏区外水域的水面风速，风速函数f(u)：

ii.对于光伏区内水域的水面风速u_pv计算公式：

u_pv＝r_s·u

其中：r_s＝u₀/u_c0；

u₀：已建成的类似光伏阵列内地面或水面风速；u_c0：已建成的类似光伏阵列外地面或水面风速；

光伏区内水域蒸发量：

E_pv＝f(u_pv)·(e_s-e_a)；

子步骤6.3.3，计算蒸发潜热：蒸发潜热L_v计算公式：

L_v＝2.5×10⁶-2.3×10³·T_s0；

子步骤6.3.4，蒸发潜热通量Q_E通过单位时间蒸发量得到：

i.光伏区外水域：Q_E＝L_vρ₀E

ii.光伏区内水域：Q_Epv＝L_vρ₀E_pv

其中：ρ₀为水体密度(kg/m³)

6.4对流热通量的计算：

子步骤6.4.1：计算波文比b：

b＝0.62P_a/1000(T_s－T_a)/(e_s－e_a)

其中：P_a为大气压力；

子步骤6.4.2：根据波文比计算对流热通量：

i.光伏区外水域：Q_H＝b·Q_E；

ii.光伏区内水域：Q_Hpv＝b·Q_Epv；

6.5水体向外长波辐射热通量E_W：水体向外长波辐射热通量计算公式：

E_w＝εσT_s ⁴；

6.6计算t时刻自由表面总热通量F：

i.光伏区外水域：

F＝I(1-α_s)+E_a-E_w-Q_E-Q_H；

ii.光伏区内水域：

F_pv＝I_pv(1-α_s)+E_apv-E_w-Q_Epv-Q_Hpv

其中，α_s为太阳短波辐射反射比；

步骤7，计算水温分布：求解物质运输方程，获得t时刻水温分布：

求解物质输运方程：

式中：c为物质浓度；D_H为水平扩散系数；D_v为垂向扩散系数；λ_d为物质衰变系数；S为源项；

步骤8，判断是否达到采集边界：判断是否T＝T_n，如果是，则结束计算，如果否则回到步骤5计算下一时段的水温分布。

2.根据权利要求1所述的方法，其特征在于，子步骤6.1.1所述的环境散射辐射量的计算方式包括如下子步骤：

子步骤6.1.1.1，计算大气层外辐射量I₀：

其中：E_sc为太阳常数；δ为赤纬角；

为当地纬度；ω为时角；λ为日-地距离变化引起大气层上界太阳辐照度的修正值：

其中：n_d为一年中的日序号；

子步骤6.1.1.2，计算晴空指数：

子步骤6.1.1.3，计算散射辐射与总辐射之比：下垫面散射辐射I_d与总辐射I之比的计算公式：

子步骤6.1.1.4，计算散射辐射：

I_d＝I·r_d。

Images