与通过拦河筑坝或引水的方式利用水头差发电的传统水轮机[1-2]相比，直接安装在河流或海洋中利用水流动能发电的水轮机不仅结构更简单，且对生态环境的破坏更小[3-4]．按转子相对于主流速度的方向可以分为水平轴[5-6]和垂直轴[7-8]两种．与水平轴水轮机相比，由于垂直轴水轮机的旋转轴垂直于水流方向，工作叶尖速比较低，具有来流适应性强，噪音污染小等优点，因此更适于安装在浅水河流．根据叶片受力方式的不同，将垂直轴水轮机分为升力型、阻力型和混合型[9-10]．Savonius水轮机[11]是阻力型垂直轴水轮机的典型代表，通过流体在迎流的凹凸叶面上形成的阻力差旋转做功产生扭矩[12]，具有工作转速低、启动力矩大、结构简单、制造及维护成本低等优势[13]．本研究设计了一种布置在浅水河岸的新型阻力型垂直轴水轮机，结构形式与Savonius水轮机类似，但部分嵌入在河岸侧壁空腔结构中．这种水轮机专门为中小型浅水河流设计，安装在河流的中层，可以有效适应流量、泥沙等多种环境因素的变化，并且由于其在低叶尖速比下工作，对水生生物的伤害也更小．通常采用数值模拟或模型试验的手段对水轮机的性能及尾流特性进行研究．KANG等[5]对水平轴水轮机进行了数值模拟，研究表明水轮机轮毂涡与叶尖剪切层之间的相互作用导致尾迹的径向曲流运动．CHAWDHARY等[6]通过数值模拟和模型试验对由三个水平轴水轮机组成的阵列的尾迹特性进行研究．OURO等[8]在水槽中对垂直轴水轮机进行模型试验，研究发现雷诺平均方程中的对流项对尾流恢复的影响比湍流输运项更大．本研究采用数值模拟法研究上述浅水河岸水轮机的性能和尾流特性，并通过模型试验加以验证．基于Fortran语言开发了一套基于有限差分法和水平集浸入边界法[14-16](level-set immersed boundary method，LSIBM)，并适用于任意复杂几何运动边界的计算流体动力学(CFD)求解器，对该浅水河岸水轮机在明渠中的流动进行了大涡模拟(large eddy simulation，LES)．通过LSIBM来描述半圆柱形空腔和水轮机的几何形状．采用耦合水平集和流体体积(coupled level-set and volume of fluid，CLSVOF)法[16]对自由液面进行精确追踪．计算了不同叶尖速比下的功率系数，并进行模型试验来验证仿真结果．研究了该水轮机的时均及瞬态尾流特性，为进一步推广该水轮机在中小型浅水河流中的应用奠定了基础．1 数值方法1.1　控制方程采用无量纲的不可压缩的连续性方程和纳维-斯托克斯动量方程模拟空气和水的流动，∇u=0；(1)∂u∂t+u⋅∇u=1ρ(-∇p+∇⋅(μS)+ρg-∇⋅τ)+f，（2）式中：u为速度矢量；t为物理时间；ρ为流体的密度；p为压强；μ为流体的动力黏度；g为重力加速度；S为应变率张量；τ为采用动态Smagorinsky模型的亚格子应力张量；f为浸入边界的附加力，以满足无滑移边界条件．1.2　水平集浸入边界法在流体和固体交界面附近的网格点上采用了LSIBM[14-15]来模拟流体与固体之间的相互作用．用水平集(LS)函数计算每个网格点到浸入边界的距离，以捕捉固体界面．LS函数ϕ在流体域为正，在固体域为负，绝对值表示网格点与浸入边界上最近点间的距离．先通过LS函数值将流场中的网格点分为流体点、固体点和IB点．再利用LS函数将IB点沿其表面法线方向投影到实体表面上，然后通过3个最近的流体点和投影点上的速度值插值得到IB点的速度，使其满足无滑移边界条件．河岸的半圆柱型空腔结构属于规则的简单几何模型，可以直接给出LS函数的表达式．但对于水轮机这种复杂几何结构，很难直接给出LS函数的计算表达式．因此提出了一种适用于任意复杂几何运动边界的固体界面捕捉方法，算法具体步骤如下：a．如图1所示将固体边界离散为STL格式的非结构化三角型网格；b．计算空间网格点到各三角形网格的距离，其中距离最小值即为该网格点到界面的距离，即为LS函数的绝对值；c．采用射线法[15]通过判断网格点与固体边界的位置关系来确定LS函数的符号，结合步骤b得到复杂结构的LS函数；d．求解类似本研究中三维旋转的水轮机的运动边界，可以根据边界的运动信息通过插值得到每一时刻的LS函数，无须在每个时间步进行重复的LS函数计算，节省了大量的计算时间．10.13245/j.hust.220324.F001图1水轮机模型示意图1.3　自由液面追踪采用CLSVOF方法[16]来捕捉空气-水两相界面．该界面捕捉算法同时解决了水平集法守恒性差及流体体积法界面精度不够高的问题，并结合两者的优点，可以对气液两相界面进行高精度的求解．密度和黏度是根据网格中水的体积分数(VOF)函数F来确定，具体为ρ=ρa(1-F)+ρwF；(3)μ=μa(1-F)+μwF，(4)式中下标a和w分别表示空气和水，对应的F分别为0和1．部分充水的网格单元的F值在0~1之间．为了捕捉两相界面，通过求解对流方程对每个时间步上的LS函数ϕ和VOF函数F进行更新，在每一个时间步利用LS函数和分段线性界面构造的算法对界面进行重构，然后根据重构的界面确定F．对ϕ和F进行更新之后，用F来校正ϕ，以保证各流体相的质量守恒．气液两相界面追踪的数值算法详见文献[16]．1.4　数值格式采用有限差分法对控制方程进行离散．空间离散基于交错网格．时间推进采用二阶龙格库塔(RK2)法．将投影法与RK2法的每一个步骤相结合以满足无散度条件．时间步长采用了基于CFL数计算的动态时间步长．利用科学计算工具包(PETSc)求解泊松方程及信息传递接口(MPI)库进行并行计算．2 数值模拟及模型试验对比分析2.1　模型试验设置在美国明尼苏达大学圣安东尼瀑布实验室进行了该水轮机的模型试验．槽道长14.6 m，宽0.9 m，水深H=0.3 m，来流速度U=0.44 m/s．如图2(a)所示，该模型为一个阻力型垂直轴三叶片水轮机，部分嵌入在右侧河岸．水轮机转子半径R=0.15 m，相应的叶片宽度B=0.1 m，高度Hb=0.1 m，转子的中心位于水深中部(y=0.15 m)．水轮机叶片与水流的最大的接触长度L=0.11 m(图2(c))．水轮机模型布置在槽道入口下游7 m处，安装在右侧河岸的一个半圆柱形的空腔中．2.2　数值模拟设置如图2所示，计算域流向(x)长度为12R，纵向(y)深度为3.6R，展向(z)宽度为3.6R．LES中的流动边界条件参考试验设置．来流速度U=0.44 m/s，自由液面高度H=0.3 m，基于U和H的槽道流雷诺数为1.32×105．由于叶尖涡和空腔壁面相互作用会产生不同程度的剪切和能量耗散，叶尖和空腔壁面的间隙大小会影响水轮机的性能．而试验中的间隙太小，仅为0.005 m，对叶尖附近的网格精度要求太高．为了平衡计算资源和网格精度，采用了3种不同的间隙值和相应的网格方案(如表1所示)．在叶片附近区域采用局部加密的均匀网格，加密区域大小如表1所示．10.13245/j.hust.220324.F002图2大涡模拟中水轮机位置及计算域示意图10.13245/j.hust.220324.T001表1模拟中计算域网格数量网格参数方案1方案2方案3间隙值/m0.015 00.010 00.007 5网格总量/1071.663.156.47叶片附近局部加密区域2.2R×1.2R×2.2R2.13R×1.13R×2.13R2.1R×1.1R×2.1R叶片附近最小网格尺度2.3×10-2R1.1×10-2R5.0×10-3R2.3　结果对比试验中控制槽道流量恒定，采用编码器对旋转轴施加不同的转速进行连续监测，得到叶尖速比(λ)的范围．再通过伺服电机控制水轮机旋转角速度(ω)，采用扭矩传感器测量对应的扭矩(T)．在模拟中，将水轮机转速设置为与试验相同的范围，并通过计算网格点上的压力和黏性力外推到水轮机表面来求得流体对水轮机表面施加的力和力矩．叶尖速比和功率系数定义为：λ=ωR/U；(5)CP=2Tω//(ρU3BHb)．(6)数值模拟中3种网格方案下的功率系数与试验结果的对比如图3所示，方案1的功率系数曲线与试验的曲线呈现相同的变化趋势，最佳的叶尖速度比约为0.3．随着网格分辨率的提高，模拟和试验之间的相对误差从方案1的16.9%减少到方案3的4.8%，因此采用网格方案3研究水轮机尾迹特性．10.13245/j.hust.220324.F003图3试验和不同网格方案下数值模拟的功率系数曲线3 尾流特性分析3.1　时均尾流特性图4(a)、(b)和(c)分别为水平横截面(左边的横截面位于水轮机中心高度y=R，右边的横截面位于y=R/3处)上平均流向速度，雷诺应力(u'w')和湍动能(k)的云图．雷诺应力表明了湍流的不稳定特性．湍动能是湍流强度的度量，具体为k=u'2+v'2+w'2/2，(7)式中：u'，v'和w'分别为对x，y和z方向的瞬时速度进行雷诺分解的脉动速度分量；⋅为时间平均符号．从图4(a)可以看出：在水轮机下游沿着水流流向，水轮机转子引起的流向速度亏损逐渐恢复．水轮机后的速度减小的区域扩展到x/R5~6的范围内，表明若布置水轮机阵列，则水轮机之间的有效间距可以略大于这个距离．从图4(b)可以看出：由于水轮机下游尾流区的涡脱落与水轮机旋转形成的周期性波动相互作用，雷诺应力在水轮机的近尾流区靠近壁面处达到峰值，即表明在这一区域湍流不稳定性最大．如图4(c)所示，当z/R=1时，水轮机尾迹产生了大量的湍动能，这与尾流区域的剪切应力的边界层相对应．近尾流区湍动能和雷诺应力的增加是由水轮机叶片周围的周期性动量通量引起的，在3.2节会进一步讨论．10.13245/j.hust.220324.F004图4不同水平横截面上流场分布云图3.2　瞬时尾流特性3.2.1　流向速度的预乘谱分析为研究水轮机尾流的周期性流动及其与时均尾流结构的关系，对流向速度进行了预乘谱分析．首先对流向速度的自相关函数进行傅里叶变换得到功率谱密度函数，乘以对应频率f得到预乘谱(G)．预乘谱用于检测引起大量湍流动能产生的主频，并沿x方向追踪湍动能能量的衰减．图5给出了流向速度的预乘谱，图中fb为叶片通过频率．可以看出：预乘谱对应的峰值在f/fb=1，这表明水轮机在尾流中产生周期性的波动流，其波动频率与fb相同．如图5所示，在流向x/L9范围内，叶片通过频率有显著的能量峰值．OURO等[8]对Gorlov型(属于升力型)的三叶片垂直轴水轮机尾流区的流向速度进行功率谱分析发现，叶片通过频率沿流向在1.5D(近尾流区)内出现了峰值，其中D为转子的直径，而叶片周围的动态失速涡引发了这个峰值．阻力型与升力型的这一差异表明，本研究中叶片产生的涡并不是造成流向速度能量峰值的唯一原因．由此可见：与升力型水轮机不同，阻力型水轮机有更显著的周期性叶片阻塞效应．阻塞效应会使来流动量发生偏转，导致暴露在水流中的叶片出现速度亏损．这也使得水轮机尾流区产生了频率为fb的周期性的波动流动，它会影响尾迹的非定常．10.13245/j.hust.220324.F005图5流向速度的预乘谱3.2.2　雷诺应力三重分解为进一步分析波动流对尾流湍流的影响，对流向速度信号进行三重分解．从速度信号u(t)=u+u˜(t)+u'(t)中分离出时间平均速度u，波动速度u˜(t)及脉动速度u'(t)分量，三重分解中的各项具体为u=1Tm∫u(t)dt；u(t)p=1N∑i=0Nu(t+iτ)    (0≤t≤τ)；u˜(t)=u(t)-u；u'(t)=u(t)-u-u˜(t)，式中：Tm为总时间；⋅p为相平均符号；τ为波动的周期；N=Tm/τ为周期数．通过比较基于三重分解(uT'wT')和雷诺分解(u'w')的侧壁沿流向的雷诺应力分布来研究周期性波动流动特性对侧壁雷诺应力的影响．图6给出在水轮机叶片中心高度处离侧壁最近的网格点的雷诺应力纵向分布．雷诺分解过程中瞬时速度与时均速度的差值由三重分解过程中的脉动和波动两部分组成(u'=uT'+u˜T)．因此，基于雷诺分解的雷诺应力在三重分解中以四项分量表示为u'w'=u˜Tw˜T+uT'wT'+u˜TwT'+uT'w˜T．从图6可以看出：u˜Tw˜T项是导致u'w'增加的最主要因素．这表明波动流动增加了近侧壁流动的不稳定性，会导致侧壁剪应力较大．而相关项u˜TwT'和uT'w˜T相互抵消，因此对u'w'的影响可忽略不计．在x/L=9后，uT'wT'和u'w'十分接近，表明波动特性消失，而脉动特性开始占主导作用，这也就解释了图5(a)的流向速度谱在x/L=9范围内检测到功率谱峰值的原因．10.13245/j.hust.220324.F006图6基于三重分解(虚线)和雷诺分解(实线)的侧壁附近的雷诺应力分布4 结论水轮机的作用引起了下游近尾流区的流向速度亏损，随着流向速度亏损逐渐恢复，速度减小的区域在x/R5~6范围内．若布置水轮机阵列，阵列之间的有效间距可以略大于这个距离．近尾流区出现了雷诺应力和湍动能增大的区域，这与通过水轮机叶片周围周期性的动量通量有关．对尾流区各监测点的流向速度信号进行了预乘谱分析，发现叶片通过频率fb处的峰值沿流向延伸到x/L=9，即尾流区的不稳定性一直延伸到x/L=9，且在尾流区检测到一种与fb相同的波动流．对尾流区靠近侧壁处各监测点的雷诺应力进行了三重分解，结果表明：叶片周围的涡脱落与水轮机产生的周期性波动流动相互作用，形成了尾迹高雷诺剪切应力区，造成了尾流区的不稳定．10.13245/j.hust.220324.F007