水平轴潮流能水轮机阵列对区域潮流场影响研究❋

袁 鹏， 刘小栋， 司先才， 王树杰， 谭俊哲， 郑志爽

(1.中国海洋大学工程学院，山东 青岛 266100； 2.青岛市海洋可再生能源重点实验室，山东 青岛 266100)

随着全球变暖等一系列由于化石燃料燃烧导致的环境问题以来，人们对可再生能源的开发和利用已经愈加重视。其中，潮流能作为一种开发潜力巨大的可再生能源，在近些年已经逐步引起了人们的关注。英国[1]、美国[2]、法国[3]等国家都已开始了潮流能的研究，在2009年欧洲海洋能中心(EMEC)颁布了第一份潮流能开发标准体系[4]并建立了较大规模的潮流能设备海试试验基地。到目前为止，各类潮流能水轮机也在不断地研究和测试中，其发展水平已接近产业化[5]。

相比于国外而言，我国潮流能技术也同样发展迅速[6]。近些年，在国内的许多水域进行了较多潮流能资源评估的研究，并对潮流能转化装置进行了研制[7]。在开发潮流能资源时，潮流能转化装置的安装可能会对其周围的自然水流模式产生影响，从而改变沉积物的分布和输运，进而影响物理环境，尤其是大型商用阵列的布放可能逐渐改变海洋物理环境，因此在进行大规模潮流能开发时有必要对潮流能转化装置阵列对水域潮流特性的影响进行研究。

以青岛斋堂岛及附近水域为研究对象，该水域的潮汐类型属于正规半日潮，最大潮差达到4.6 m，平均差约为2.8 m。其中，在斋堂岛东南部水域水深为35～40 m且水底相对较为平坦，大潮期间峰值流速能超过1.9 m/s，是我国北方比较有代表性的潮流能资源区。近些年来，许多关于潮流能的研究都是围绕该水域进行的：在2012年，李华军等[8]通过Flux方法对其琅琊台海峡的潮流能资源进行了评估；同年，邵萌等[9]对在该水域建立500 kW海洋能独立电力系统示范工程进行了规划；在2013年，史廉博等[10]基于FVCOM建立三维水力学模型并对水域潮流特性进行分析；在2014年，纪合盼等[11]通过实测水文数据对潮流能发电站的建设进行初步规划；在2015年，陈娅玲等[12]采用动量损失的方法模拟水轮机，并对该水域的水轮机阵列方式进行了研究；在2017年，林杰等[13]基于Delft 3D研究了该水域的漩涡等因素对潮流能提取可能带来的影响。从现有文献来看，目前的研究多集中于对装置发电能力的优化以及区域潮流能资源的评估，而对潮流能装置对海洋环境的影响研究还很少。

本文基于Delft 3D-Flow建立斋堂岛及附近水域的三维水力学模型，并通过实测数据对模型进行验证。经过验证后的模型，在合理确定该水域中水轮机阵列放置位置[14]的基础上，利用等效[15]方式模拟水平轴潮流能水轮机阵列，对加入水轮机阵列后的水力学模型进行计算，分析水轮机阵列对区域流场的水位和流速等影响，为未来在此水域开发潮流能资源提供参考。

1 潮流场水力学模型

1.1 建立模型

本次研究采用Delft 3D的Flow模块建立水力学模型。基于浅水方程和Boussinesq假设，该模型采用有限差分法求解斜压下的N-S输运方程，既可以用于二维模型，也可以用于三维模型计算[16]。其主要基于三组控制方程，分别是：

质量守恒方程：

(1)

动量守恒方程：

(2)

(3)

输运方程：

(4)

其中:ζ表示水位(总水深)；d表示到参考平面的当前水深(净水深)；U、V分别代表在x，y方向上的速度分量；Q代表单位面积上质量源强度；f为科氏力参数；vh是动态水平涡流粘度；ρ0是参考密度；ρ′是不规则密度；τsx、τsy是作用在海平面上的风压分量；τbx、τby为底部的剪应力分量。在输运方程中，c代表盐度或者温度；Dh是水平涡流消散度；λd代表一阶衰减过程；g为重力加速度；R为单位面积上的源项。

本次研究的模拟水域为青岛斋堂岛及其附近水域。斋堂岛隶属于山东省青岛市，位于琅琊台东南方向海中，地理坐标为北纬35°38′，东经119°55′(见图1)。

为提高模型输出结果的计算精度和减少计算机资源利用，在建模过程采用网格嵌套(Nesting)方法，分别建立2个水域范围和网格分辨率不同的模型(见图2，图中黑色区域为大范围水域网格范围，红色区域为小范围网格范围，蓝线为陆线)。

图1 斋堂岛地理位置

大范围水域所建立的模型由其开边界上通过OTIS提供的8个潮汐调和常数(M2、2、N2、K1、K2、O1、P1、Q1)进行驱动，考虑到OTIS的输出的潮汐调和常数的分辨率(1/30(°))，需要建立足够大小的水域计算网格范围保证足够数量的调和常数作为边界条件，从而确保模型仿真结果的可靠性。网格分辨率由水平方向上的水深数据分辨率决定，网格尺寸选取600 m×600 m并在斋堂岛区域局部加密到120 m×120 m。

小范围水域模型的建模水域范围由本次研究所需的研究水域范围决定，其驱动使用大范围水域模型输出的水位时间序列作为输入。网格尺寸选取36 m×36 m并局部加密到18 m×36 m以保证后续对水轮机阵列等效模拟的精度。

图2 模型网格

计算时大范围水域采用二维水力学模型，小范围水域采用三维水力学模型。其中，小范围的水力学模型在垂直方向上将水流分为10 sigma层(每层厚度为10%的水深)并采用k-ε湍流模型进行计算。在确保满足Courant条件的前提下，大小模型的时间步分别设置为1和0.1 min。利用所建立的水力学模型模拟一个半月潮的潮汐周期并且不考虑风、浪等影响，模拟时间选择为17 d，其中前2 d作为模型由静水状态开始运行的启动时间。

1.2 潮流场水力学模型验证

为确保计算结果的有效性，分别通过实测的水位数据对大范围水域模型进行验证和实测流速流向数据对小范围水域模型进行验证。

1.2.1 大范围模型验证 在大范围水域模型中，选取日照港、董家口港和斋堂岛潮位站3个站点2013—2014年的部分实测水位数据进行验证。其中大范围模型的验证点坐标以及部分验证结果分别见表1和图3。

通过实测水位数据与模型模拟结果进行对比：日照港水位验证的平均水位误差为5.58%，局部最大误差为13.06%；斋堂岛的平均水位误差为8.00%，局部最大误差为15.10%；董家口港平均水位误差为9.94%，局部最大误差为13.86%。

表1 潮位验证点经纬度

1.2.2 小范围模型验证 在大范围水域模型经验证后，将其输出的水位时间序列作为小范围水域模型的边界条件进行输入， 然后将小范围水域模型的输出结果与实测流速流向数据进行对比，并在此基础上调节小范围水域模型的参数从而提高模型输出结果的精确度。该水域:4个流速验证点的坐标和位置分别见表2、图4。

图3 斋堂岛潮位站验证结果

流速验证Velocity verification站点1Site 1站点3Site 3站点4Site 4站点5Site 5经度Longitude/°E119°55′39.9″119°55′46.5″119°55′39.2″119°56′8.79″纬度Latitude/°N35°37′27.8″35°37′27.3″35°37′18.5″35°36′12.88″

根据Rahman A等的研究[17]，水力学模型的输出流速受水底摩擦系数影响较大。而本次建立水力学模型的水域，因为缺乏斋堂岛附近水域的精确水底摩擦数据，此次仿真采取统一的水底摩擦系数。在此，分别选用Manning=0.017、0.019和0.021对该区域流速验证进行对比。图5为小范围水域模型的流速误差验证。

将模型的仿真结果和实测流速进行对比可知(见图5)：1号站点流速验证离散型较大且整体仿真流速偏大，可能由于其离岛较近受到岛屿尾流所引起的漩涡影响较大[13]且水深变化剧烈，所以仿真精度较低；4号站点和5号站点处水深变化较为平缓且距离陆地相对较远，流速验证准确性较高；3号站点由于仍处于斋堂岛尾流的影响范围之内，相对于4号和5号站点流速验证准确性略低。

此外，根据图5对比水底摩擦Manning系数不同时模型输出各个站点的均方根误差(RMSE)和相关系数(r)的数值，由此可以看出可知当Manning取0.019时，各站点的RMSE相对较小且r较大，因此流速验证精度最高。因此，在后续的仿真模拟中，模型的Manning系数都选用0.019。

图4 流速验证点位置

图5 各站点流速验证

2 阵列对区域潮流特性改变研究

2.1 水轮机阵列放置位置选择

潮流能水轮机阵列的位置通过TSE方法[14]来确定。TSE方法主要用于水深有限的水域中对水轮机位置进行选址，该方法很好地兼顾了区域潮流能能量强度和水深因素，并引入了潮流能可开发度(Tidal Stream Exploitability, TSE)因子来表征某区域中适合放置水轮机的程度，TSE因子的求解公式为：

(5)

式中：V0和h0分别为该水域的典型流速和水深，在这里分别取1 m/s和20 m；Vf和Ve分别为该水域涨落潮期间的流速；h为涨落潮期间的局部水深；ξ为补偿函数，其计算公式为：

(6)

(7)

(8)

其中：h为涨落潮期间的水深；Δh为最大潮差；h1和h2分别为水轮机安装水深补偿的最大上下限，在这里分别取2、5 m。

根据本次研究所建立的小范围水域模型可得到该水域涨落潮阶段的流速分布(见图6)。仿真结果显示，涨落潮阶段在斋堂岛东南部区域和其与大陆间的琅琊台海峡处流速较大。在大潮期间，涨潮流速能超过2.0 m/s，落潮流速能超过1.8 m/s。

图6 大潮涨落潮阶段流速分布

采用TSE方法可计算出该水域的TSE因子分布(见图7)。在斋堂岛东南部区域，TSE因子较大，比较适合放置规模较大的水轮机阵列；而在大陆与斋堂岛间的海峡处，虽然涨落潮期间流速较大，但水深相对较小，不适合放置较大规模的水轮机阵列。因此，考虑到实际水深等因素，将水轮机阵列初步选定布放于斋堂岛东南部水域，水轮机阵列位置如图7黑线内所示，其中阵列放置区域的面积为540 m×360 m。

(黑线内部为水轮机阵列的放置区域。Black line is the placement area of the turbine array.)

图7 斋堂岛附近水域TSE因子分布图
Fig.7 TSE index around waters of Zhaitang Island

2.2 水轮机阵列设置方式

在本次研究中，水轮机阵列中的水轮机都假设为水平轴潮流能水轮机。在流场模型中加入潮流能水轮机依研究的目的不同有不同的处理方法。其中，一些研究通过修改水力学模型中计算网格局部网格节点上的水底摩擦系数来等效水轮机阵列所带来的水动力影响[18]，不过此方法多用于二维水力学模型[19]；另外一些研究通过在水力学模型中计算网格的局部节点上添加动量损失方程来仿真潮流能水轮机对水域产生的影响[12]；此外，也有一部分研究通过修改水力学模型中的局部湍流模型系数来模拟水轮机对潮流场带来的改变[20]。

本次研究通过上述第二种方法，采用动量损失的方法。该方法的具体过程如下：

在水流通过潮流能水轮机时，会产生一定程度的动量损失，根据公式(2)和(3)，其动量方程分别可变为：

(9)

(10)

其中，Mx和My分别为水力学模型在U和V两个方向上的动量源项。

在Delft 3D中，可以通过UDF文件多孔盘(porousplate)来实现对其水力学模型中局部网格节点上动量方程的源项修改，以实现对潮流场中水轮机的模拟。该方法已通过实验与Delft 3D中的仿真结果进行了一定验证[15]并在近年的不少研究中得到了应用[19]。

图9 多孔盘的在计算网格中布置方向

使用多孔盘仿真潮流能水轮机时，除了需要计算出水流通过潮流能水轮机时产生能量损耗的大小，还需要考虑在实际情况下的水轮机迎流方向与建立的水力学模型中计算网格方向不匹配的问题，因此需要在水轮机放置位置上的计算网格中的U和V两个方向上同时设置多孔盘(见图9)。多孔盘中参数设置的计算公式为：

(11)

其中，动量耗散参数Closs设置的计算公式，在U方向上为：

(12)

(13)

Au=∑At|cosθ| 。

(14)

同样，在V方向上：

(15)

(16)

Av=∑At|sinθ| 。

(17)

其中：γu和γv分别代表在U和V方向上单位网格上的动量损耗比例；CT为水轮机的轴向力系数，本次仿真假设CT为定值并且选取CT=0.8[21]；Au和Av分别代表在U和V方向上单位网格上的水轮机叶片扫掠面积；At代表在单位网格上的水轮机的叶片上的扫掠面积；Δx和Δy分别代表单位网格节点上在U和V方向上的距离；Δz代表每sigma层的高度；n代表多孔盘所在的层数。

2.3 水轮机阵列对区域水位流速影响

本次研究，所选择布放水轮机阵列的区域平均水深在35 m左右，根据EMEC的标准[4]，水平轴水轮机叶片最底端需要高出水底5 m以上，水轮机叶片最顶端需要低于该水域天文潮水位最低点5 m以上，因此选用叶片直径为18 m的水轮机，其中水轮机转子中心位置在距离海床高度15 m处。

在所选区域放置4排水轮机，每排间距离为10D(D为水轮机转子直径)；每列中相邻水轮机转子中心间距离分别选用1D、2D、3D、4D和5D，采用交错布置的方式。在水轮机阵列布放面积确定的条件下，以上各阵列布置方式中可布放水轮机台数分别为80、40、28、20和16台。分别将未放置水轮机和放置不同排布方式的水轮机阵列的潮流场模拟结果进行对比可得到的半月潮大潮期间水位和流速矢量对比如图10所示。

仿真结果显示：涨落潮期间，水位会在水轮机阵列迎流方向上有一定程度的上升，且每经过一排水轮机后会有一定的下降，该现象随阵列中水轮机密度的增加而愈加显著。其中，当一定面积阵列中水轮机密度达到最大时，其对区域水域的改变最大不超过3 cm，且随其与阵列距离的增加影响幅度逐渐减小。另外，阵列对区域水位影响范围较大，将会影响到岛屿与大陆间海峡处的水位(见图10(a))。

相比对水位的影响，水轮机阵列的布放密度对流速影响十分明显，其所带来的最大流速下降在尾流影响区域能够超过0.45 m/s，且影响距离超过2 km(见图10(a))。而当选择较为稀疏的水轮机阵列时(行间距大于4D时)，其对潮流场的影响较小，流速改变不到0.1 m/s(见图10(e))。此外，在该水域加入水轮机阵列后将会对由岛屿尾流所引起的漩涡[13]带来一定改变，此种改变所带来的影响需要在未来进一步研究和评估。

((a) 列间距为1D,水轮机台数共80台；(b) 列间距为2D,水轮机台数共40台；(c) 列间距为3D,水轮机台数共28台；(d) 列间距为4D,水轮机台数共20台(e) 列间距为5D,水轮机台数共16台。(a) 80 turbines with 1D horizontal spacing； (b) 40 turbines with 2D horizontal spacing； (c) 28 turbines with 3D horizontal spacing； (d) 20 turbines with 4D horizontal spacing； (e) 16 turbines with 5D horizontal spacing.)

图10 不同阵列布置方式下的水轮机阵列对区域水位和流速影响

Fig.10 Effect of turbine array on water level and velocity in different layout

3 结语

本文以Delft 3D的Flow模块为基础，建立斋堂岛及其附近水域的三维水力学模型。通过TSE方法确定了水轮机阵列的布放位置，并将水轮机对流场的影响等效成为动量损失的方式并组成不同排列方式的阵列加入到水力学模型中，从而对流场进行仿真并研究阵列对潮流场所产生的影响。

仿真结果显示，在斋堂岛东南部水域放置一定规模的水平轴潮流能水轮机阵列将会对该区域的流场带来一定影响：阵列会对较大范围水域的水位带来一定改变，且影响程度随阵列中水轮机密度的减少而降低，其中最大水位改变不超过3 cm；对区域流速的影响随阵列中水轮机密度的增加而愈加显著，当密度较大时阵列所引起的流速下降在其尾流区域能够达到0.45 m/s，影响距离可延伸至其下游2 km处。