一种风电场风机尾流动态耦合模拟方法

本发明属于计算风工程领域，更具体地，涉及一种风电场风机尾流动态耦合模拟方法。

背景技术：

在风电场中，风经过风机之后会产生尾流区，也就是一个风速衰减区，这种现象称为尾流效应，这会导致处于风电机组下游的风机的发电功率显著降低。因此，开展风机尾流的研究，对提高风电机组发电功率有着重要意义。随着计算机技术的发展和提高，cfd(computationalfluiddynamic，计算流体动力学)的方法广泛用于风机尾流的数值模拟，国内外许多学者的研究大多基于此展开。

早些年的研究中，闫海津等对风力机进行全模型化数值模拟，得到整个流场的速度分布，压力分布和湍流分布等结果。但是，进行全模型建模时，风机叶片贴体网格划分难度大，网格需求量高，计算效率明显降低。近年来，任会来等提出了带旋转的致动盘模型，以致动盘代替风轮圆盘，该方法基于致动盘模型和bem(bladeelementmomentummethod，叶素动量理论)，根据不同位置处的翼型特点，计算轴向和切向诱导力，然后将诱导力沿着翼展方向微小圆环内做平均，以此获得不同位置处的轴向与切向诱导力。该方法在cfd数值模拟中引入体积力源项，实现风机对风场的扰动效果，如此既可以避免全模型划分网格造成的计算效率降低，又可以同时考虑风机所受轴向和切向诱导力，从而该方法的数值模拟精度具有一定的可靠性。

但是，该方法的局限性在于，在cfd计算流场的迭代过程中，致动盘所在的网格单元添加的体积力源项(即风机阻力项)始终为一定值。而实际上，风机基于bem理论划分的微小圆环对应体积力的大小与致动盘前的来流风速有关，处于风电机组下游的风机来流风速会受上游风机尾流的影响，在计算迭代过程中不断变化，所以各圆环对应体积力源项也应因此而变化。因而，亟需一种耦合计算该流场和体积力源项的方法，实现二者的联动，修正数值模拟风机对该流场的扰动作用，更加真实的计算风机尾流流场的分布特性。

技术实现要素：

针对现有技术的以上缺陷或改进需求，本发明提供了一种风电场风机尾流动态耦合模拟方法，其目的在于，基于cfd数值模拟，通过在计算域的指定位置处提取致动盘各圆环所对应的单元网格，然后对应添加阻力项，并根据流场每一步迭代计算得到的各圆环处对应的来流风速，计算并更新对应单元网格的阻力项大小；该方法能够耦合计算计算域的流场和风机对该流场的扰动效果，更加真实的反应风机的尾流分布情况。

为实现上述目的，本发明提供了一种风电场风机尾流动态耦合模拟方法，包括如下步骤：

s1、将计算域划分为多个网格单元，进而建立计算域网格模型，同时获取风机相关参数；

s2、将风机致动盘划分为多个径向等距的圆环，根据风速分别计算各圆环对应的轴向力和切向力，进而将轴向力和切向力对应分配到圆环包含的每个网格单元，得到各网格单元对应的阻力项，并依此建立风速与阻力项的对应关系模型；

s3、基于计算域网格模型，对风机尾流进行数值模拟，模拟时将阻力项添加到对应的网格单元，并根据变化的风速和风速与阻力项的对应关系模型动态计算更新阻力项，完成风机尾流的动态耦合模拟。

作为进一步优选的，轴向力fn和切向力ft根据下式计算得到：

其中，ρ为空气密度，n为风机叶片数量，dr为圆环的径向宽度，c为圆环处对应的风机叶片弦长，cn和ct分别为轴向推力系数和切向阻力系数，ω为叶片上的合速度。

作为进一步优选的，叶片上的合速度ω根据下式计算得到：

其中，vn为圆环包含网格单元的平均轴向风速，vt为圆环包含网格单元的平均切向风速，ωr为圆环处绕轮毂中心旋转的线速度，vt+ωr为圆环的相对切向风速。

作为进一步优选的，轴向推力系数cn和切向阻力系数ct计算步骤如下：

(1)根据圆环的平均轴向风速和相对切向风速计算入流角

(2)根据下式，由入流角得到轴向推力系数cn和切向阻力系数ct：

其中，先由入流角和桨距角β得到攻角α：进而根据攻角α得到圆环对应的升力系数cl和阻力系数cd。

作为进一步优选的，所述s1中，先根据风电场所在地形区域确定计算域，对计算域进行网格划分得到计算域网格；然后获取计算域内的高程数据，并对高程数据进行插值处理；最后将插值处理后的高程数据与计算域网格结合，得到计算域网格模型。

作为进一步优选的，对计算域进行非结构化网格划分，划分后每层网格数量相等且为三角形网格单元，每个网格单元整体呈不规则三棱柱状。

作为进一步优选的，将轴向力和切向力对应分配到每个网格单元时，首先对圆环包含的单元总体积作平均，得到单位体积的轴向力和切向力，然后将单位体积的轴向力作为网格单元x方向的阻力项，将单位体积的切向力进行y、z方向分解后，分别作为该网格单元y、z方向的阻力项，最终得到每个网格单元对应的x、y、z方向上的单位体积的阻力项。

总体而言，通过本发明所构思的以上技术方案与现有技术相比，主要具备以下的技术优点：

1.本发明基于cfd数值模拟，对致动盘方法进行了改进，通过在计算域的指定位置处提取致动盘各圆环所对应的单元网格，然后对应添加阻力项，并根据流场每一步迭代计算得到的各圆环处对应的来流风速，计算并更新对应单元网格的阻力项大小，实现了在计算风电场的流场过程中动态添加风机阻力项；本发明与流场耦合计算的方法，能够更加准确的模拟风机与流场的相互作用，并且更加真实的反应风机的尾流作用效果。

2.本发明给出了根据风速计算轴向力和切向力，并进一步得到阻力项的具体计算方式，且该动态添加风机阻力项的方法避开了使用冗杂的建模方法来模拟风机对风场的作用，显著提高了数值模拟风电场的效率和效果，此外可根据需求自定义调整风机的类型、数量、位置、轮毂高度以及风机致动盘划分圆环的数量等参数信息，具有操作简单、灵活便捷，功能全面、实用性强等优点。

3.本发明对计算域进行了非结构化网格划分，且划分的网格单元整体呈不规则三棱柱状，对应关心区域有着更高的网格精度，并且非结构化网格具有生成速度快、数据结构简单、网格质量好等优势，容易实现边界的区域拟合，相比结构化网格更适用于cfd计算。

附图说明

图1为本发明实施例风电场风机尾流动态耦合模拟方法流程图；

图2为本发明实施例中两台风机轮毂高度附近xy截面上阻力项分布云图；

图3a和图3b分别是本发明实施例中上游风机和下游风机致动盘yz截面阻力项分布云图；

图4为本发明实施例中两台风机数值模拟所得的流场的速度分布云图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。

本发明实施例提供的一种风电场风机尾流动态耦合模拟系统，如图1所示，包括如下步骤：

s1、将计算域划分为多个网格单元，进而建立计算域网格模型，同时获取模型和风机相关参数；

s2、将风机致动盘划分为多个径向等距的圆环，分别计算各圆环对应的轴向力fn和切向力ft，进而将轴向力fn和切向力ft对应分配到圆环包含的每个网格单元，得到各网格单元对应的阻力项，并依此建立风速与阻力项的对应关系模型；

s3、基于计算域网格模型，对风机尾流进行cfd数值模拟，即通过数值求解流体运动的控制方程，得出流场在连续区域上的离散分布，从而近似模拟风机尾流分布情况；模拟时将阻力项添加到对应的网格单元，并在每次迭代中，根据风速变化和风速与阻力项的对应关系模型动态计算更新控制方程的阻力项，以修正数值模拟风机对流场的扰动，从而实现阻力项与流场耦合计算，直到达到计算迭代收敛，完成风机尾流动态耦合模拟。

进一步的，步骤s1具体包括以下子步骤：

1.1选取研究地形区域，确定计算域，确定计算域中心的经纬度坐标，从地理空间数据云gis系统中导出包括该计算域的高程数据文件，该dem数据文件有30m精度的分辨率，选取计算域时，确保计算域能够包含足够多的数据点以反映该计算域内的地形环境信息。

1.2对计算域进行非结构化网格划分，每层网格数量相等且为三角形网格单元，划分的网格单元整体呈不规则三棱柱状，对应关心区域有着更高的网格精度。

1.3将dem高程数据文件导入服务器端gis平台，用自编程序对相邻的数据做一定的插值处理，然后结合已划分好的计算域网格，在计算域网格底端创建地形，计算底面网格端点被纵向抬升或降低的高度，最终得到可用于计算的复杂地形cfd模型，即计算域网格模型。

1.4将风机相关参数信息储存在独立文件，便于在程序运行过程中读取并使用这些参数，然后用于致动盘体积力源项的计算；具体的，风机相关参数包括：多台风机在计算域中的位置坐标以及该风机的轮毂高度等可修改的预输入风机信息；风机的翼型参数信息，具体包括将致动盘划分得到的同心圆环的半径r，以及该半径处对应的翼型弦长和桨距角；叶片气动文件，具体包括风机叶片沿翼展方向上的攻角和对应的升力系数和阻力系数。

进一步的，步骤s2具体包括以下子步骤：

2.1将风机致动盘划分为多个径向等距的圆环，确定并提取出各圆环所包括的网格单元，各圆环的径向宽度为dr；根据叶素动量理论，作用在风轮平面dr范围的叶素上的轴向力fn和切向力ft分别为：

其中，ρ为空气密度，n为风机叶片数量，c为dr处所对应的风机叶片弦长，cn和ct分别为轴向推力系数和切向阻力系数，ω为叶片上的合速度。

2.2叶片上的合速度ω由以下公式计算得到：

其中，vn为圆环处的平均轴向风速，vt为圆环处的平均切向风速，ωr为圆环处绕轮毂中心旋转的线速度，vt和ωr之和为圆环处的相对切向风速。

具体的，在进行cfd数值模拟迭代计算时，每次迭代都提取致动盘网格单元对应的xyz三个方向的风速(以地形高度方向为z轴建立空间直角坐标系)；对于由致动盘划分的每个圆环，提取圆环包含单元的x方向的风速后作平均，以此作为该圆环的轴向风速，即vn；提取并合成单元yz方向的风速并求和作平均，即vt，然后加上该圆环绕轴旋转的线速度，即ωr，得到该圆环的相对切向风速，即vt+ωr，轴向风速和相对切向风速二者合成即可得到合风速ω。

2.3计算出圆环对应的轴向风速和相对切向风速后，由下式计算入流角

其中，先由计算得到攻角α，β为桨距角，对于不同的风机，不同叶素处对应的桨距角已给定，在该程序中由已储存的叶片翼型参数文件查得；进而根据攻角α由叶片气动文件插值得到圆环对应的升力系数cl和阻力系数cd，由以下公式合成轴向的推力系数cn和切向的阻力系数ct：

2.4基于以上步骤便能在每一迭代步中计算得到致动盘各圆环对应的轴向力fn和切向力ft；将这两种力具体对应分配到每一个网格单元上时，首先对圆环包含单元总体积作平均，得到单位体积的轴向力和切向力，然后再将相对于致动盘的切向力根据每个单元相对于圆环圆心的位置关系进行分解，最终得到每个网格单元对应的xyz方向上的单位体积的阻力项。

以下为具体实施例：

1、选取研究地形区域，确定计算域，根据gis高程数据建立复杂地形下的计算域网格模型；

1.1选取经度116.2°、纬度39.8°为目标计算域的底面中心原点坐标，选定长宽高都为800m的计算域模型，并从gis系统中获取该计算域的地形高程dem数据文件；

1.2使用网格划分程序对该计算域进行网格划分，每层网格数量相等且为三角形网格单元，每个单元整体呈不规则三棱柱状；

1.3将dem数据文件导入到服务器端，然后使用自编程序对已经划分好的计算域网格模型的底层网格根据dem数据文件进行z坐标的抬升或降低，插值模拟出复杂地形的cfd模型，为了后续数值模拟计算的稳定与收敛，在计算前将模型缩尺100倍。

2、将风机的相关参数信息储存在独立文件，在程序运行过程中读取并存储风机信息，然后用于致动盘体积力源项的计算；

2.1将需要布置的两台风机的位置、轮毂高度单独存储一个文件，程序读取后用于后续计算；

2.2分别将计算域相关的参数、风机的翼型相关参数、风机的叶片气动参数单独存储一个文件，程序读取后用于后续计算。

3、计算风机轮毂中心在计算域中的实际坐标，并以此为基点确定风机致动盘各圆环包括的网格单元；

3.1对于该复杂地形的cfd模型，边界条件设置如下：入口采用速度入口，设置为10m/s，垂直于速度入口，方向同x轴正方向，出口为压力出口，侧面采用对称边界条件，底面地形采用壁面；

3.2风机轮毂中心在计算域中位置即为致动盘的圆心坐标，而风机轮毂中心实际z坐标为风机轮毂高度和该点地形高度之和，计算出各台风机所在点的地形高度，加上对应轮毂高度可得其在计算域中的实际z坐标；

3.3以各风机轮毂中心为基点，沿风机径向将致动盘划分为八个圆环，然后遍历整个计算域的所有单元网格，判断单元的中心是否在对应的圆环内，从而确定并提取出各圆环所包括的单元网格。

4、根据流场迭代过程动态计算更新致动盘包含单元网格对应的轴向和切向诱导力，从而获得该风机与流场耦合计算后的数值模拟结果，如图2所示；

4.1提取各风机各圆环处的x方向平均风速作为该圆环的轴向风速，提取各风机各圆环处的yz方向合速度再平均后再加上对应的线速度作为该圆环的相对切向风速，进而计算出入流角以及合速度ω；

4.2入流角对应减去文件存储的桨距角β可得攻角α，设定空气密度ρ＝1.293kg/m³，风机叶片数量n＝3。根据叶片气动文件数据插值得到升力系数cl和阻力系数cd，然后进行分解与合成得轴向升力系数cn和切向阻力系数ct；

4.3计算得到圆环所受轴向和切向的诱导力，然后对圆环包含单元总体积作平均算得单位体积的轴向力和切向力，将该轴向力作为单元x方向的阻力项，如图3a和3b所示，将切向力进行yz方向分解后再分别作为该单元yz方向的阻力项，最后，添加该阻力项到cfd数值模拟计算过程，并且每一迭代步都更新该阻力项的值；

4.4最终计算收敛后，数值模拟计算这两台风机的尾流分布情况，如图4所示。

本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。