在沙漠、戈壁和荒漠地区，干旱少雨的环境正是太阳能富集区的天气特征之一，地理和气候因素决定了二者分布的一致性．修建光伏电站需要大量的土地资源，沙戈荒地区逐渐成为修建大型集中式光伏电站的主要区域，但沙漠光伏电站的生态效应却少有人关注．光伏厂区内与地表切应力有关的地表侵蚀和堆积，与光伏板周围流场有关的光伏板表面覆尘等，对光伏电站的安全性和发电量有着极大影响．研究[1-4]表明：覆尘将最大电流显著降低6.9%～16.4%，由积尘的沉积热辐射损失可达到15%．光伏阵列第一排与第二排光伏板更容易积尘，Ali Salari等[5]通过数值模拟得到光伏组件表面的灰尘沉积密度从0 g/m2增加到8 g/m2，电效率降低了26.36%；赵伟平等[6]通过数值研究得到对于粒径和年平均风速分别小于30 µm和7 m/s的区域，适当增加倾角可以减少灰尘积聚．对于风载荷的研究，Abiola等[7]通过风洞试验得到了光伏板上下板面的压力系数分布；Jubayer，You等[8-10]通过数值模拟与风洞试验分析了单排和多排光伏阵列各排光伏板的升力、阻力及详细的流场特征．黄斌等[11]通过风洞试验，建立了沙尘暴气候下的风沙流场，结果表明：在风沙流时光伏组件的风载荷特征及较低风速下，含沙量随高度呈指数级下降，当风速增加时，含沙量呈反S型下降．华亚平等[12]对光伏电站内外遥感图像中检索到的气温、地表温度和反照率数据的气象观测数据进行分析，指出光伏电站内的气温、地表温度和反照率均低于电站外，且随季节表现为冬季明显，夏季不明显[13]．切应力可以理解为动量通量，光伏阵列可以视为粗糙元，气流经过光伏阵列时光伏板可以分配一部分动量通量，作用在地表的动量通量即为地表剪切应力总量肯定会减少[14]．τs/τ(τs为存在粗糙元的地表应力，τ为没有粗糙元时的地表切应力)是粗糙元动量分配理论中的重要参数．亢力强等[15]通过风洞试验得到了不同植株作为单一粗糙元时的动量分配规律．可以将光伏阵列理解为粗糙元并承担了一部分的动量通量，导致τs在后方区域降低，从而起到固沙的作用．地表切应力还与输沙量息息相关，地表剪切应力与输沙量成正比关系．光伏电站内局部区域产生的风蚀积现象与地表切应力的不均匀分布有关，即存在气流加速区域增加沙尘颗粒的运动，而气流减速区减少沙尘颗粒的运动，既沉积区．因此，在保证沙漠光伏电站发电效率的同时应让光伏阵列有一定的防风固沙能力．目前对于沙漠光伏电站固沙效应的研究较少，大多通过测量光伏阵列对风速的削弱，得到光伏阵列具有固沙能力的结论．野外对于光伏电站固沙效应定量的研究，测量难度较大，风洞试验须对光伏板缩比，但对沙尘颗粒却无法应用相似理论进行缩比，所以试验结果可能会与实际结果存在出入．本研究开展光伏电站内风场数值模拟，主要考察与风蚀积最相关的地表切应力的大小和分布，并利用沙尘颗粒的流体起动公式给出了不同风向下的保护区范围，分析了光伏板周围的流场特征，为进一步分析沙尘的运动即覆尘问题提供基础理论．1 数值模拟方法1.1　计算模型如图1所示，研究对象为西北腾格里沙漠南部边缘地区某光伏电站内成排的光伏阵列，为了减少计算量，取每排光伏板宽4 m、长18 m，即宽度和实际尺寸一致，长度缩减．10.13245/j.hust.250631.F001图1西北腾格里沙漠沙漠边缘地区某光伏电站与简化模型为了考察不同排数对地面固沙效应并减少计算量，考虑到个别光伏电站由于地形的限制，光伏板在计算域而对于多排光伏板简化至三排光伏阵列．第一排代表光伏电站内前端光伏阵列，第三排代表末端光伏阵列，第二排代表中间光伏阵列．光伏板前沿据地面0.35 m，后沿距地面2.74 m，光伏板与地面倾角α=37°，θ为风向角[16]．如图2所示，当光伏板前沿迎风时风向角为正，当后沿迎风时方向角为负．通过旋转整个光伏阵列模拟不同风向下的工况，图2(a)中虚线表示旋转45°的光伏阵列为例．X方向为流向，光伏板前沿和后沿之间的距离为行道，其中第一排与第二排之间为第一行道，第二排与第三排之间为第二行道．10.13245/j.hust.250631.F002图2光伏阵列位置前沿、后沿、风向角1.2　网格划分及计算设置图3为计算域及光伏板网格图，采用非结构网格对整个流场进行网格划分．非结构网格具有良好的结构适应性．10.13245/j.hust.250631.F003图3计算域及光伏板网格所选用的湍流模型是Realizable k-ε，壁面模型选用标准壁面函数模型．为保证计算的合理性保证y+在30~600的范围之内可通过公式y+=ρμτywall/μ(ρ为流体密度；μτ为流体壁面摩擦速度；ywall为第一层网格高度；μ为运动黏度)计算第一层网格高度，保证网格划分的合理性．整个计算域为一个150 m×70 m×20 m的长方体，通过旋转光伏阵列模拟不同风向来流，在地面添加了17层近地面网格，第一层网格高度为2 mm．使用Fluent2020R2行计算，边界条件设置如下：计算域入口为速度入口，出口为压力出口，两侧及顶部设置为敞口对称边界．底部和光伏板均为无滑移壁面．速流场入口初始速度为对数廓线μ(z)=μ*κlnzz0，式中：κ为冯卡门系数；z0为床面粗糙度；μ*取0.3 m/s．使用SIMPLE算法，离散格式选用二阶迎风格式，稳态计算，估算可得地表切应力应为0.11 Pa．1.3　保护区分析方法沙粒的起动分为流体起动和冲击起动，流体起动作为沙粒起动过程中的第一阶段，是判断是否起沙的重要依据，但是流体起动的机制较为复杂，是近地表湍流、微观粗糙度和地表切应力共同作用的结果[17]．沙粒被吹离地面返回沙床撞击其他沙粒会造成冲击起动，击溅更多的沙粒在空气的裹挟下跃移运动，连锁反应进而形成风沙流，所以控制沙粒流体起动是固沙的关键．临界地表切应力是判断沙粒是否会发生流体起动的准则．沙粒起动的临界摩阻风速为u*t0(ds)=a1ρp-ρaρagds+b1ρads，换算得到临界起沙的地表切应力τt0=ρaμ*t02，(1)式中：ρa为空气密度，取值1.225 kg/m3；ρp为沙粒密度，取值2 650 kg/m3；ds为沙粒直径；在土壤干燥情况下a1=0.11；b1=3×10-4 kg/s2；g为重力加速度．以腾格里沙漠的分布最多的120 μm沙粒为标准[19]，得到了临界起沙的摩阻风速．根据式(1)可以计算出临界起沙的地表剪切应力为0.051 6 Pa，小于临界地表切应力则认为地表的沙尘颗粒不发生流体起动．这些区域将它定义为保护区．2 计算验证与Abiola-Ogedengbe等[7]的风洞实验和Jubayer等[8,10]的数值模拟结果进行对比，设定与Jubaye一致的条件(光伏板长7.29 m，宽2.48 m，倾角25°，风向角90°)，光伏板上下板面的压力系数(Cp)分布如图4(a)所示．如图4(b)所示(τs和τ分别为地表切应力实际值和理论值)，设置了三种网格节点50×108，100×214，200×428，发现100×214兼顾计算效率与计算准确性．10.13245/j.hust.250631.F004图4压力系数及网格无关性验证3 计算结果及分析3.1　光伏阵列流场图5展示了θ=±90°时，光伏阵列y=0 m截面的X方向的速度分布云图．图6为θ=±90º，±45°时光伏阵列在z=0截面速度云图．10.13245/j.hust.250631.F005图5θ=±90º时光伏阵列y=0截面X方向速度云图(色标单位：m/s)10.13245/j.hust.250631.F006图6θ=±90º，±45°时光伏阵列在z=0截面速度云图(色标单位：m/s)θ=90°时，光伏阵列板下和板间内风速明显降低，在第一排光伏板上方出现绕流加速现象，而下方产生较大的回流涡，近地风速由于狭管作用会出现局部增大，主要集中在第一排光伏板的板下区域和第一行道．在光伏阵列最后一排板后也会形成较大的涡，但未靠近地表．在光伏阵列内流线总是沿流线相反方向，第三排光伏板后沿后方区域也会形成较大的涡．θ=-90°时，同样地，在光伏阵列内风速削弱明显，在第一排光伏板前后出现复杂的流体涡系，板下由于阻塞效应出现较小的回流涡，前沿后方会出现旋转方向相反的涡对，靠近地表由于狭管效应流体加速且涡强较大，其上方的涡范围较大但涡强较小；这与Jubayer等[7]的数值模拟得到流场云图一致．θ=±90º时，气流由于绕流会在光伏阵列行道产生涡，风驱动沙粒会使得光伏阵列行道产生风蚀积现象．但在光伏阵列内的风速却远低于周围环境风速．结合图5和图6综合分析可知：θ=±90º时，光伏场内呈现三维的涡系结构，主要表现为前端的两排和靠近两侧的区域；而θ=45º时风速减少的区域集中在上风向的一侧；θ=-45º时风速减少的范围较大．3.2　地表切应力图7为θ=±90°，±45°时沙漠地表切应力τs分布等值线云图．10.13245/j.hust.250631.F007图7θ=±90º，±45°时光伏阵列地表剪切应力分布等值线云图(色标单位：Pa)可以发现：θ=90°时第一排光伏支架处最大，第一排板下区域及第一行道存在较大τs分布，在光伏阵列两侧也存在较大的切应力分布．郝璐等[11]通过数值模拟得到：θ=90°时，第一排至第五排光伏板的沙尘最大沉积率分别为18.89%，12.35%，9.62%，6.83%和5.71%，这与地表切应力分布具有一致性．所以θ=90°时，沙漠光伏电站防沙的重点区域为第一排光伏板下与第一排和第二排行道．在光伏阵列后方τs存在明显降低．θ=-90°时，在光伏板前沿下方区域地表切应力较大，由于光伏板的阻塞作用，支架附近没有出现τs增大的现象．θ=-90°时整个光伏阵列板前沿下方地表均会出现剪切应力局部增大，在光伏阵列两侧会产生气流加速，导致τs增大．所以在主风向θ=-90°的沙漠光伏电站，应更注重于对光伏阵列板前沿下方地表风蚀的防治，光伏阵列第一排与最后一排是防沙的主要区域．由图7(c)可以看出：τs最大区域为光伏板支架区域，在光伏阵列迎风一侧地表存在较大的切应力分布．在后方部分区域依然存在切应力减小的区域，但是较大值的切应力分布面积较小，在考察的四组风向中，此风向应注重迎风侧后沿位置的局部防沙．θ=-45°时，切应力最大的区域在光伏板前沿下方地表，而不是在支架附近，这与θ=45°时不同．θ=45°时地表切应力增大的区域明显大于θ=-45°，此风向下应注重光伏板两侧的防沙治沙．由以上分析可知：不同风向下地表切应力分布差异较大，应根据光伏电站当地的主风向采取相应的防沙治沙措施．为了考察沿风向光伏电站不同典型位置地表切应力的变化，取y=0.0，4.5，9.1 m位置处的地表剪切应力分布绘于图8(a)～(c)，x的取值范围为-40~75 m，X=115 m，虚线表示光伏板前沿和后沿的位置．可以看出，地表剪切应力在第一排光伏板前沿地表达到最大．除了前沿位置处，在光伏阵列前后方及光伏阵列内部区域，三个不同位置的τs/τ均小于1．这说明光伏阵列的铺设可降低地表的切应力．10.13245/j.hust.250631.F008图8θ= ±90º，±45°时光伏阵列不同位置处的剪切应力分布由图8(a)可以看出：在y=0 m处，θ=90°时τs/τ在光伏阵列前逐渐降低，在光伏阵列前方降至最低，然后在第一排光伏板前沿位置显著增大，τs/τ可达2.8后又显著降低至0.1，在第二排与第三排光伏板下方τs/τ会有小幅增大，但其值小于1．θ=-90°时，在第一排光伏板前沿下方增大，τs/τ增大为3．在第一行道内急剧降低，降低至0.07，在第二排光伏板后均小于0.5，在第三排光伏板前沿下方地表增大，增大至1，降低后该位置地表切应力逐渐恢复．由图8(b)可知：在y=4.5 m处，θ=90°地表切应力分布在第一排光伏板前沿下方地表急剧增大，τs/τ增大为2.8．在第三排前沿下方地表切应力增大，τs/τ增大为1.2．风向角θ=-90°时，地表切应力在第一排光伏板前沿下方增大，τs/τ增大为3.2，在第一行道内急剧降低至为0.13，在第二排前沿处会轻微增大，但小于1．在第三排光伏板前沿下方地表切力增大，增大至1.8，降低后该位置地表切应力逐渐恢复．由图8(c)可以看出：在y=9.1 m处，θ=90°时第一排光伏板前沿下方地表切应力增大，τs/τ增大为2.4左右，第三排的增大量总是大于第二排．θ=-90°时地表切应力分布在第一排光伏板前沿下方地表切应力急剧增大，τs/τ增大至3.4后在第一行道内急剧减小，τs/τ减小至0.35，在第二排前沿下方地表切应力增大，τs/τ增大至1，后又在第二行道内τs/τ减小为0.34．在第三排光伏板前沿下方地表切应力增大，τs/τ增大至2.3，在后方τs/τ减小为0.6后，地表切应力逐渐恢复．由以上分析可知，光伏电站内地表切应力的变化呈现二维分布．在图8(d)中，y0=0 m位置处，θ=45°时的地表切应力分布高于θ=-45°，这表明45°来流时光伏阵列内地表更容易发生“风蚀”．由图8(c)和(d)可以看出：θ=45°时，切应力局部增大的面积要大于-45°；θ=±45°时，在流场中的该位置处的切应力增大的区域总是在光伏阵列板前沿下方地表．与θ=±90°不同，θ=±45°时在该位置，光伏阵列每一排前沿地表切应力均会增加．θ=45°时增加量最大的是第三排光伏板，τs/τ可达2.9，第一排增幅最小，τs/τ可达1.8，第二排可达2.6．θ=-45°时，在该位置第一排前沿地表切应力增幅最大，τs/τ可达2，其次是第三排τs/τ可达1.8，第二排τs/τ可达1.75．θ=±90º，±45°时第一至第三排前沿地表切应力分布见图9，由图9(a)可知：θ=90°时第一排光伏板前沿下方存在较大的切应力分布，τs/τ可达到2．但第二排和第三排光伏板除了光伏板边缘处τs/τ大于1，其他区域τs/τ均小于1．在y/Y=0.2~0.8区间，第三排板前沿下方高于第二排，而两侧区域却是第二排其值高于第三排．由图9(b)可知：三排光伏阵列板前沿下方地表切应力分布存在相似的规律，即沿着板沿至中间位置切应力分布减至最小值，两侧近乎对称分布．第一排和第三排前沿地表存在较大的切应力分布，均高于代表中间行列的第二排，第一排τs/τ均在2.5以上，第三排均在1.2以上．这说明在主风向为θ=-90°时，沙漠光伏电站第一排与第三排前沿下方以及两侧是防沙的重点区域．由图9(c)可知：θ=45°时光伏阵列各排板前沿切应力顺着前沿地表切应力先增大后逐渐减小，其中第三排的切应力分布高于第二排和第一排，所以主风向为10.13245/j.hust.250631.F009图9θ=±90º，±45°时第一至第三排前沿地表切应力分布θ=45°的光伏电站光伏阵列前沿地表均存在比±90º时较大的切应力分布，须在此区域做防沙处理．由图9(d)可知：θ=-45°时光伏阵列板前沿地表切应力分布比θ=45°复杂，在迎风侧存在先增大后减小的趋势，但在另一侧存在先增大后减小非对称分布的趋势，迎风侧第三排的分布曲线高于第一排，第一排高于第二排，在另一侧第一排高于第三排高于第二排．这说明第二排的前沿地表切应力总体较低．整个光伏阵列前沿地表τs/τ均高于1.5，即表现为电站内整体增大的情况．3.3　保护区分析对于运行一段时间已经出现风蚀积现象的沙漠型光伏电站，由于风场的筛选作用，电站内出现不同的粒度分布[19]．研究颗粒流体起动，地表切应力可作为判断颗粒是否发生流体起动的依据．选取沙尘粒径为腾格里沙漠区的标准粒径正态分布，其平均值为120 μm[20]，定义低于临界地表切应力的区域为保护区．虽然目前起沙机制还存在争议，但沙粒的临界流体起动摩阻风速可以一定程度反映沙尘颗粒开始运动的动力学行为．图10为θ= ±90º，±45°三排光伏阵列形成的保护区(虚线为光伏阵列的大体位置)．10.13245/j.hust.250631.F010图10θ= ±90º，±45°三排光伏阵列形成的保护区θ=90°时，保护区较θ=-90°宽度小，但长度却更长．θ=±45°保护区的位置明显减小，θ=45°时保护区位于光伏阵列的行道内，θ=-45°时在光伏阵列前方和行道内会形成保护区．建议在对沙漠光伏电站进行防沙设计过程中可以在图10虚线内空白区域采取一定的措施．4 结论a．θ=90°时，地表切应力最大的区域为第一排光伏板的支架附近，其次为第一排光伏板前沿下方，增大180%．θ=-90°时地表剪切应力最大的区域为第一排光伏板前沿下方，最大可增大250%．光伏阵列前方与后方区域地表切应力τs/τ均小于1．b．θ=±45°时，平行于光伏阵列后方为地表切应力减小的区域．θ=45°时地表切应力最大的区域出现支架附近，θ=-45°地表切应力最大的区域为光伏板前沿迎风一侧下方．θ=45°风向y0=0位置处，在光伏板前沿下方地表切应力的值和增大区域的范围均大于θ=-45°风向．c．θ=±90º时，在光伏板前沿下方地表和光伏板两侧地表切应力最大，θ=90°时增大100%，θ=-90°时增大200%．θ=±45°时在迎风侧光伏板前沿下方地表区域切应力最大：θ=45°时，增大200%，θ=-45°时增大180%．d．在θ=90°风向下，第二排与第三排前沿下方地表切应力会降低，且τs/τ均小于1．θ=-90°时光伏阵列每一排前沿下方的切应力分布曲线趋势相同．在光伏板的两侧切应力最大，越接近光伏板中部，地表切应力越小．θ=±45°时，光伏板迎风一侧存在较大地表切应力．e．光伏阵列有一定的固沙能力，在光伏阵列前后方会形成保护区，保护区的面积和位置会随风向发生改变．可以根据不同风向来流的保护区形状对沙漠光伏电站做防沙设计．10.13245/j.hust.250631.F011