CN111651896B - 一种基于实际风速及实际复杂地形的风电场流场计算方法 - Google Patents

Abstract

本发明涉及风电场流场计算技术领域。技术方案是：一种基于实际风速及实际复杂地形的风电场流场计算方法，包括如下步骤：步骤1，利用Python语言实现基于相关性分析的自动插补测风数据以及风资源分析；步骤2，利用快速傅里叶变换得到基于实际风速的动态模拟风速函数；步骤3，利用DEM高程数据得到风电场真实地形三维结构图；步骤4，根据步骤2‑3所得的结果，利用LBM与LES相结合的新型CFD数值计算方法，采用动态自适应方法进行格子排布与分裂，对不同区域采用不同的空间离散尺度，并结合浸入边界法对动态旋转边界进行仿真，实现整个风电场空气动力场的仿真。该方法可靠性高，有助于提高风电场微观选址的准确性。

Description

一种基于实际风速及实际复杂地形的风电场流场计算方法

技术领域

本发明涉及风电场流场计算技术领域，具体涉及一种基于实际风速及实际复杂地形的风电场流场计算方法。

背景技术

风能作为一种清洁并且储量丰富的可再生能源在全球范围内的开发利用已经得到了很大发展，其中风力发电是现代社会中与人们生产生活息息相关的一种风能利用的主要手段。作为风电项目开发的第一步，科学有效的风资源评估对提高风电场经济效益有着不可忽视的作用，风电场的微观选址对于风电场项目效益也颇具影响。

受天气变化、仪器故障等客观条件影响，测风设备记录的观测数据可能有部分失真甚至缺失等现象发生。因此在进行风能资源评估之前，首先要对拟选场址实际测风数据进行预处理，完成相关的检验。

计算流体力学(Computational Fluid Dynamics，CFD)方法通过数值计算求解控制方程进而研究风力机空气动力学特性，从其计算结果可以获得旋转风力机诱导流场的详细信息，现在已经成为风机气动特性计算、尾流研究和风电场布局优化等问题的主要研究手段之一。现有针对风电场流场的研究大多基于传统CFD方法，对于边界入流条件的设定一般选用均匀入流，精细度不够，传统CFD计算方法对复杂地形和风力机动态旋转边界处理困难。此外，地形模型与实际地形的误差也是影响数值计算精度的重要因素。

数值模拟方法中新兴的格子玻尔兹曼方法(Lattice Boltzmann Method，LBM)是CFD领域的一个强有力分支，近年来受到国内外许多学者的重视并取得了开拓性进展，相较于一些传统数值模拟方法具有诸多优势，如编程简单、易于处理复杂边界特别是运动边界的问题等。

与传统方法不同，格子玻尔兹曼方法的基础是分子动力学理论，但是并不关注每个分子的运动行为，而是通过分布函数表示包含大量分子的流体粒子运动特征。格子玻尔兹曼方法与大涡模拟(Large Eddy Simulation，LES)的结合具有良好的适用性，在可及计算资源下获得高雷诺数湍流流场并可用于求解非稳态运动边界的尾流结构，此外，其中的格子自适应排列方法还可针对复杂几何边界和尾流结构进行自动优化。

发明内容

本发明所要解决的技术问题是针对上述背景技术的不足，提供一种基于实际风速及实际复杂地形的风电场流场计算方法；该方法应具有可靠性高、有助于提高风电场微观选址的准确性的特点，从而极大提高风电场的经济效益。

本发明为实现上述发明目的采用如下技术方案：

一种基于实际风速及实际复杂地形的风电场流场计算方法，包括如下步骤：

步骤1，利用Python语言实现基于相关性分析的自动插补测风数据以及风资源分析；

步骤2，利用快速傅里叶变换并基于余弦函数拟合得到基于实际风速的动态模拟风速的时域函数；

步骤3，利用DEM高程数据得到风电场真实地形三维结构图，并获取风力机位置，构建考虑真实地形和风力机的流场计算域；

步骤4，根据步骤2-3所得的结果，利用LBM与LES相结合的新型CFD数值计算方法，采用动态自适应方法结合浸入边界法进行格子排布与分裂，对不同区域采用不同的空间离散尺度，并结合浸入边界法对动态旋转边界进行仿真，实现整个风电场空气动力场的仿真。

进一步的，步骤1中的相关性分析，采用最小二乘法降低损失函数，自动插补的序列包括风速和风向数据。

进一步的，步骤2中的快速傅里叶变换，对秒级测风数据进行离散傅里叶变换分析，以频域分解各分量的幅值作为参考依据，提取频谱中的主要分量，拟合出时域上的函数作为模拟风速序列。

离散傅里叶变换的具体方法如下所述：

对连续信号x(t)u(t)以T_s为抽样时间间隔进行均匀理想抽样，得到有限长的离散信号即因果序列x(n)，对其求得的傅里叶变换称为序列的离散时间傅里叶变换(Discrete-Time Fourier Transform，DTFT)，记为X(e^jθ)，其中θ＝ωT_s：

离散时间非周期序列的傅里叶逆变换为：

将式(1)对应的频谱进行离散化，取一个主值区间，在周期2π内划分为N等分，则

由式(1)和式(2)可得离散傅里叶变换对：

快速傅里叶变换(Fast Fourier Transform，FFT)算法在DFT基础上进行改进，当N为合数时有N＝n₁n₂，则原长度为N的DFT可以分解为两个长度分别为n₁和n₂的DFT，以此类推可以继续进行分解成更多长度的DFT来计算。

进一步的，在步骤3中，DEM数据通过Matlab编程对其进行插值加密并绘制地形三维图像，地形数据以三维文件输出。

进一步的，在步骤4中，LBM方法基于Boltzmann-BGK方程的离散形式而建立，包括空间、时间和速度的离散，流场在空间上被划分为格子区域，时间上按一定的时间步长进行离散，其速度被简化为有限维数的速度空间。

进一步的，在步骤4中，所述LES采用的亚格子湍流模型为WallAdapting LocalEddy(WALE)模型。

进一步的，在步骤4中，所述自适应方法结合浸入边界法用较大尺度的格子进行模拟，并对特定位置进行局部细化。

与现有技术相比，本发明所达到的有益效果是：本发明结合大涡模拟、自适应格子排列以及浸入边界法，发展了基于LBM-LES方法的风电场空气动力场多尺度仿真模型；通过Python语言编程高效准确地实现测风数据插补和风资源分析，基于快速傅立叶变换进行频谱分析并实现了脉动风速的建模，避免均匀来流条件对风电场流场模拟造成的误差；通过Matlab编程建立风电场复杂地形的三维实体模型，并将风力机三维结构图置于相应位置构建流场计算域。基于LBM-LES数值计算方法进行流场计算，为风电场中风机排布提供支撑。本发明预测的可靠性高，实现了实测风速和实际复杂地形共同作用下的风电场流场仿真，可得到较真实的风电场空气动力场及风力机尾流时空演变规律，有助于提高风电场微观选址的准确性和精细化程度，极大提高风电场的经济效益。

附图说明

图1为测风数据插补程序实现流程图。

图2为实测风速序列与拟合曲线对比图。

图3为实际地形三维建模技术路线。

图4为基于实际地形数据生成三维地形结构的示意图。

图5为自适应网格加密技术示意图。

具体实施方式

下面结合附图对本发明的优选实施例对本发明的技术方案进行详细说明：

(一)基于Python语言及线性回归分析的测风数据插补；

受天气变化、仪器故障等客观条件影响，测风设备记录的观测数据可能有部分失真甚至缺失等现象发生。因此在进行风能资源评估之前，首先要对拟选场址实际测风数据进行预处理，完成相关的检验。

原始数据为实施例风电场不同高度测风塔测风数据，格式为*.xls文件，包含日期时间、风速大小及风向等信息。根据图1所示测风数据插补程序实现流程图，利用Python编写相应程序，读取格式为*.xls文件的测风数据，建立参考测风数据和目标测风数据两个数据表，生成两个数据结构类型变量，变量包含日期时间、风速和风向三个序列，分别标记出变量中风速、风向两列数据包含的空值；提取出两个测风塔风速数据中非空的值，用剔除空值风速序列的构造出两个列表类型变量。

读取数据模块的输出为4个变量：参考测风数据和目标测风数据的数据结构变量、清洗过的参考风速序列和目标风速序列的列表变量。回归分析模块的输入是为参考测风塔和目标测风塔两个风速数据的列表型变量，构造出回归对象，用最小二乘法找出最佳的线性模型拟合两个数据集之间的线性关系，绘制数据的散点图和拟合直线，求出常数项、回归系数、判定系数并得到回归方程，用以上参数构造一个字典变量作为函数返回值。

数据插补模块的输入参数有4个：目标测风数据和参考测风数据的数据结构变量以及线性模型的常数项和回归系数，每个数据结构变量的时间序列是完整的，以时间列的元素个数作为完整数据量。分别遍历风速值列、风向值列，识别出目标测风塔数据中之前被标记的空值，根据索引查找并判断参考数据同期数据是否有效，如果有效便可对相应的目标塔数据插补，其中风速值用线性模型插补，风向值用参考塔同期数据直接替换；如果无效则对相应的目标数据去除标记还原空值。插补完成后生成一个新的测风数据数据结构变量，同样有日期时间列、风速数据列和风向数据列，统计处理后风速数据和风向数据的完整率并输出，将新生成的数据结构转换为*.xlsx文件格式输出。

(二)基于快速傅里叶变换的风速序列拟合；

在对原始测风数据进行插补得到完整风速数据后，本发明对秒级风速序列进行离散傅里叶变换分析，分析频域计算的结果，提取对幅值贡献较大的余弦分量，构成一条新的风速曲线，拟合成与测得风速曲线较为接近的关于时间的函数，目的是使得风场流场的模拟结果更接近于实际情况。

离散傅里叶变换的具体方法如下所述：

对连续信号x(t)u(t)以T_s为抽样时间间隔进行均匀理想抽样，得到有限长的离散信号即因果序列x(n)，对其求得的傅里叶变换称为序列的离散时间傅里叶变换(Discrete-Time Fourier Transform，DTFT)，记为X(e^jθ)，其中θ＝ωT_s：

离散时间非周期序列的傅里叶逆变换为：

将式(1)对应的频谱进行离散化，取一个主值区间，在周期2π内划分为N等分，则

由式(1)和式(2)可得离散傅里叶变换对：

快速傅里叶变换(Fast Fourier Transform，FFT)算法在DFT基础上进行改进，当N为合数时有N＝n₁n₂，则原长度为N的DFT可以分解为两个长度分别为n₁和n₂的DFT，以此类推可以继续进行分解成更多长度的DFT来计算。

在本实施例中，流场数值计算时长设定为十分钟，选取一段数据量为N的序列进行拟合，N的大小为600。首先需要对序列进行隔直处理，即将每个数值都减去所有数据的平均值，得到一个新的序列，使得零频分量的幅值变小，在后续的功率谱分析中可以减少直流能量对低中频信号谱分析的影响。通过程序读取实测风速序列，由于无周期性且均匀离散无须再进行采样，对其进行FFT运算，输出的结果是一个长度和原序列相等且元素为复数的一维数组，即为原序列的傅里叶系数。对运算结果取模并进行归一化处理，得到的序列即为频谱图的幅值。

原序列时域上抽样间隔为1s，则频谱的周期等于时域上的抽样频率f_s＝1,离散傅里叶变换将频域上每个周期划分为N等分，每根谱线的宽度为f₀＝f_s/N，所以横坐标频率为kf₀，(k＝0,1,2,……，N-1)。编程处理频谱的幅值，将归一化的模值和对应的复数形式傅里叶系数放在同一个数据结构中，由于数据结构中不便存放复数，故将实部虚部分别存放在两个序列中，并新增两列用于存放计算得到的频率和相角。频率等于频谱图的横坐标kf₀，k即为该变换系数在数据结构变量中对应的行号；相角由系数的实部虚部值根据式(4)计算得到。按照系数模值大小对数据结构降序排列，计算并填补频率列和相角列，处理后的部分数据结构如下表所示：

上表所示的数据结构变量值中，index列表示频域横轴的序号，mod列存放傅里叶系数的模值，com_r列和com_i列分别存放了傅里叶系数的实部和虚部，f列表示频率，ang列用于存放计算得到的对应相角。提取零频分量并选取模值较大的部分傅里叶系数，通过傅里叶逆变换计算得到一个时域上的函数，得到隔直处理后的函数为：

将拟合得到函数的图像与FFT变换前的实测风速序列图像进行对比，如图2所示，实线为实测风速序列，虚线为拟合所得的曲线，且能用特定的函数式表达。可以发现拟合曲线较好的反应了原风速序列的数值特征，峰值、振幅、频率以及变化趋势等信息都与实际风速曲线较吻合。

(三)复杂地形三维建模；

图3为实际地形三维建模技术路线。所要实现的是将实际地形数据转化三维模型数据(参见图4)。图4中的A图为地理空间数据云下载所得风电场实际地形的DEM高程数据图，选取框中区域，通过以下复杂地形三维建模方法，将该区域转换成Matlab能识别的基于实际地形的三维模拟图像(如图4中的B图所示)，并将其转换为三维模型作为文件输出用于数值模拟。

地理空间数据云是由中国科学院计算机网络信息中心建设和维护的一个数据及模型共享云平台，提供丰富的地学遥感数据资源，以中国区域为主，覆盖全球地理范围。

本实施例根据实际地形的DEM高程数据生成三维地形结构图的具体步骤如下：

选取风电场项目，根据经纬度从地理空间数据云平台查找到相应区域的DEM文件，选择一块指定大小的区域导出，以*.xyz的格式输出。*.xyz文件包含的是点云数据，存放了所有数据点的三维坐标，包括位置及高度，通过后续编程处理可以生成三维模型。三维坐标数据文件可以以文本形式打开，部分数据内容如下所示：

得到实际地形的点云坐标数据文件后，通过Matlab编程将地形数据转换为可移植的实际地形三维模型数据。将包含点云数据的*.xyz文件作为输入值，根据列提取xyz三个数列上的坐标数据，先对数据点进行插值，生成实际地形的三维曲面，并将三维地形表面数据写入到*.stl文件中，从而将其转换成三维模型并作为文件输出，供后续模拟使用。实现相应部分的代码如下：

(四)基于LBM-LES的风电场流场数值模拟

根据(二)中所得隔直处理后的风速函数，选取平均风速较高(6.94m/s)的实测风速序列，对风速序列插补、频谱分析并经过拟合得到的实测风速随时间的变化函数如下：

将该风速函数作为入流条件。

根据(三)中所得实际地形的三维模型，建立风场的计算域模型，并根据以上实测风速结果，基于实测风速对风场空气动力场进行仿真。选取风场计算域的尺寸为2070m×1600m×1830m，格子排列采用自适应方法，选用的湍流模型为大涡模拟(LES)地形处的网格分辨率为11.5m，计算域最大网格尺寸为46m。计算时间步长取0.1s，总计算时间为600s，使用CFD方法进行数值模拟。

图5为自适应网格加密技术示意图。在数值模拟过程中，求解精度与网格尺度的细分成正相关。若对全局网格采用同样精细的网格，会造成计算资源的极大浪费。根据解的变化和需要，对网格尺度进行自动调整，以提高数值计算效率及精度。如图5所示，在风力机叶片周围的流体域，自动加密网格以提高求解精度，而在逐渐远离叶片的流体域，网格尺度逐渐增大，这是由于较粗的网格也能满足相应区域的求解精度。该方法通过在数值模拟软件中设置完成。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进、替换或变型，都应涵盖在本发明的保护范围内。

Claims (5)

1.一种基于实际风速及实际复杂地形的风电场流场计算方法，包括如下步骤：

步骤1，利用Python语言实现基于相关性分析的自动插补测风数据以及风资源分析；

步骤2，利用快速傅里叶变换得到基于实际风速的动态模拟风速函数；

步骤3，利用DEM高程数据得到风电场真实地形三维结构图；

步骤4，根据步骤2-3所得的结果，利用LBM与LES相结合的新型CFD数值计算方法，采用动态自适应方法进行格子排布与分裂，对不同区域采用不同的空间离散尺度，并结合浸入边界法对动态旋转边界进行仿真，实现整个风电场空气动力场的仿真；

所述步骤1中的相关性分析，是采用最小二乘法降低损失函数，自动插补的序列包括风速和风向数据；

所述步骤2中的快速傅里叶变换，是对秒级测风数据进行离散傅里叶变换分析，以频域分解各分量的幅值作为参考依据，提取频谱中的主要分量，拟合出时域上的函数作为模拟风速序列；

所述步骤3中的DEM数据，是通过MATLAB编程对其进行插值加密并绘制地形三维图像，地形数据以三维文件输出。

2.根据权利要求1所述的基于实际风速及实际复杂地形的风电场流场计算方法，其特征在于：所述步骤4中的LBM基于Boltzmann-BGK方程的离散形式而建立，包括空间、时间和速度的离散，流场在空间上被划分为格子区域，时间上按时间步长进行离散，其速度被简化为有限维数的速度空间。

3.根据权利要求2所述的基于实际风速及实际复杂地形的风电场流场计算方法，其特征在于：所述步骤4中的LES采用的亚格子湍流模型为Wall Adapting Local Eddy模型。

4.根据权利要求3所述的基于实际风速及实际复杂地形的风电场流场计算方法，其特征在于：所述步骤4中的自适应方法用格子进行模拟，并对特定位置进行局部细化。

5.根据权利要求4所述的基于实际风速及实际复杂地形的风电场流场计算方法，其特征在于：所述离散傅里叶变换的具体方法如下所述：

对连续信号x(t)u(t)以T_s为抽样时间间隔进行均匀理想抽样，得到有限长的离散信号即因果序列x(n)，对其求得的傅里叶变换称为序列的离散时间傅里叶变换(Discrete-TimeFourier Transform，DTFT)，记为X(e^jθ)，其中θ＝ωT_s：

离散时间非周期序列的傅里叶逆变换为：

将式(1)对应的频谱进行离散化，取一个主值区间，在周期2π内划分为N等分，则

由式(1)和式(2)可得离散傅里叶变换对：

快速傅里叶变换(Fast Fourier Transform，FFT)算法在DFT基础上进行改进，当N为合数时有N＝n₁n₂，则原长度为N的DFT可以分解为两个长度分别为n₁和n₂的DFT，以此类推可以继续进行分解成更多长度的DFT来计算。

Images