Wind power plant wake flow prediction method based on LBM-LES

Abstract

The invention relates to the technical field of new energy wind power, in particular to a wind power plant wake flow prediction method based on LBM-LES, and the method comprises the following steps: (1) obtaining wind resources and wind power plant topographic data, and constructing a wind power plant computational domain; (2) applying boundary conditions and material parameters: applying the corresponding boundary conditions to the computational domain according to the data obtained from the step (1), and endowing the material parameters; (3) embedding an LES large eddy simulation model in the LBM through the effective viscosity in the collision model, and meanwhile, dispersing a lattice Boltzmann equation and solving the lattice Boltzmann equation; the method has the advantages that a complex macroscopic flow motion equation does not need to be solved, instead, the lattice Boltzmann equation of the mesoscale is solved, and compared with a traditional method, the solving efficiency is greatly improved; secondly, a wind turbine is simulated by adopting a rotary actuating disc method, so that a wind power plant wake flow simulation result is closer to a real wind field, and the calculation precision is remarkably improved.

Description

Technical Field

[0001] The present invention relates to the technical field of new energy wind power, and in particular to a method for predicting a wind farm wake based on LBM-LES. Background Art

[0002] Wind farm microsite selection is a technology that optimizes the layout of wind turbines to maximize the overall power generation efficiency of wind farms by accurately analyzing factors such as wind speed, wind direction, and terrain. This technology collects a large amount of meteorological data, such as wind speed, wind direction, temperature, and humidity, by building wind towers in areas with abundant wind resources. Combined with terrain data and geological conditions, it evaluates and analyzes potential wind turbine construction locations, further considers grid access conditions, road transportation conditions, etc., and ultimately determines the optimal wind farm microsite selection plan.

[0003] In the process of micro-site selection of wind farms, the wake of wind turbines and the wake interference effect are factors that must be considered. The wake effect will cause the wind speed behind the wind turbine to decrease, causing the output power of the downstream wind turbines in the wind farm to drop significantly; and the high-intensity turbulence and wind shear in the wake will also significantly increase the fatigue load of the unit and reduce the service life of the unit. The wake interference between units will further strengthen the above adverse effects. Studies have shown that the wake effect in wind farms will cause the average annual power generation of wind farms to decrease by about 15%, and in extreme cases, the power generation loss can reach 30%.

[0004] Therefore, how to achieve efficient and accurate prediction of wind farm wake effects in large-scale and complex engineering scenarios is a key technical issue that needs to be urgently solved in wind power development.

[0005] The existing wind farm wake effect prediction technology often uses analytical wake models, such as Jensen and Frandsen models, which are superimposed with pre-calculated flow field results without wind turbines. The accuracy of simulating wake effects is poor, which will greatly reduce the accuracy of wind farm micro-site selection results. Secondly, when calculating the flow state of wind farms, finite volume methods are often used. This method requires the construction of high-quality three-dimensional meshes, which is particularly difficult when dealing with complex terrains. It often requires a lot of time and manual intervention to optimize the mesh quality. In addition, this method has a high computational cost, especially when simulating large-scale wind farms or long-term series data. The amount of calculation is huge, resulting in a long simulation process, which is difficult to meet the needs of fast iteration solutions in the micro-site selection of wind farms. To this end, the present invention aims at many shortcomings of existing wind farm wake prediction technology, such as poor calculation accuracy and low efficiency, and proposes a wind farm wake prediction method based on LBM-LES, where LBM-LES is lattice Boltzmann method-large eddy simulation. Summary of the invention

[0006] In view of the above problems, the present invention provides a method for predicting the wake of a wind farm based on LBM-LES. The specific scheme is as follows: A wind farm wake prediction method based on LBM-LES includes the following steps: (1) Obtain wind resources and wind farm terrain data and construct the wind farm calculation domain; (2) Applying boundary conditions and material parameters: Based on the data obtained from step (1), apply corresponding boundary conditions to the wind farm calculation domain and assign material parameters; (3) By replacing the effective viscosity in the collision model, the LES large eddy simulation model is embedded in the LBM, and the lattice Boltzmann equation is discretized and solved; (4) Calculate the aerodynamic force of the wind turbine rotor according to the wind turbine parameters in the wind farm, and use the Gaussian projection method to map the aerodynamic force into the volume force on the cell; (5) Calculate the density and velocity of the macroscopic flow field in the calculation domain of the wind farm and obtain the wind turbine wake velocity distribution information.

[0007] The acquired wind resource data includes wind speed, wind direction, turbulent velocity, etc. At the same time, the wind farm calculation domain is constructed through the wind farm terrain data; and according to the wind resource data and wind farm data, the corresponding boundary conditions are applied to the calculation domain and the material parameters are assigned. The wind turbine parameters in the wind farm include at least the wind turbine array and model parameters; the wind turbine rotor aerodynamic force is calculated based on the acquired data.

[0008] This technical solution numerically simulates the flow state in the wind farm by solving the discrete lattice Boltzmann equation, and transmits the flow information, such as wind speed and wind direction, to the wind turbine rotating actuator disk model. At the same time, the LES large eddy simulation model is coupled in the LBM; in this process, efficient and accurate simulation of the wind turbine wake effect is achieved, and the power and thrust of the wind turbine are calculated.

[0009] Furthermore, the specific method of coupling the LES large eddy simulation model in the LBM described in step (3) is as follows: The LES large eddy simulation model is coupled to the LBM through the sub-grid eddy viscosity model. The sub-grid eddy viscosity model is formed by superimposing the sub-grid viscosity coefficient considering the small-scale eddy dissipation on the original viscosity coefficient to form an equivalent turbulent viscosity. , thereby replacing the original viscosity coefficient and changing the collision model; (1) In the formula, is the original viscosity coefficient, in Pa·s, For The subgrid viscosity coefficient considering small-scale eddy dissipation is superimposed on the basis, and the unit is Pa·s; is the equivalent turbulent viscosity, in Pa·s; C s is the Smagorinsky constant, i.e. the Smagorinsky model constant, with a value range of 0.1~0.2; is the filtering scale of the filter, taking the cube root value of the grid volume in LBM, that is, , Δ1, Δ2, Δ3 are the lengths of the unit cell in different directions, in m; defined by the following formula: (2) (3) In the formula, is the modulus of the strain rate tensor, is the strain rate tensor, in s -1 ; u i is the macroscopic fluid velocity, in m / s; i, j The value is 1, 2, or 3, corresponding to the x, y, and z directions respectively.

[0010] Furthermore, the process of discretizing the lattice Boltzmann equation in step (3) is: The continuous lattice Boltzmann equation, that is, equation (4), is discretized to obtain the discrete form of the lattice Boltzmann equation, that is, equation (5): (4) (5) In the formula, express t The velocity at position x is c i The number of microscopic particles is also called the distribution function; t is time, unit is s; c i is the velocity of microscopic particles, in m / s; F i is the force exerted by the wind turbine on the flow field, in N; ▽ is the gradient operator; As the collision operator, the Bhatnagar-Gross-Krook operator is used, that is, ; In the formula, is the local equilibrium distribution function in discrete velocity space is the weight coefficient defined in the selected discrete velocity model; is the fluid density in kg / m 3 ; u is the macroscopic fluid velocity, in m / s; c s is the lattice sound velocity, dimensionless; tis the time step, in seconds, is the relaxation time, dimensionless.

[0011] Furthermore, in step (3), the lattice Boltzmann equation is evolved and solved using the relaxation migration method to obtain the distribution function after the update time step; the specific steps are as follows: First, the collision solution step is performed: (6) In the formula, The distribution function after solving for collisions; Then perform the migration evolution solution step: (7).

[0012] Furthermore, in step (4), the method for calculating the aerodynamic force of the wind turbine rotor is as follows: According to the momentum blade element theory, the blade is discretized into N blade element segments along the span direction. Assuming that the blade element segments do not interfere with each other, each blade element segment can be regarded as a two-dimensional airfoil extension. Furthermore, based on the incoming wind speed of the blade element segment and the aerodynamic parameters of the airfoil, the aerodynamic load of each blade element segment is calculated. The aerodynamic force of the blade element segment at any radial position r of the wind rotor is expressed as: (8) In the formula, F r is the aerodynamic force of the blade element at any radial position r of the wind rotor, in N / m; dr is the blade element length, in m; C L is the airfoil lift coefficient, C D is the airfoil drag coefficient, which can be obtained by looking up the table based on the airfoil information; e L is the direction of the lift unit vector, dimensionless, e D is the direction of the drag unit vector, and its unit is dimensionless; is the air density, in kg / m 3 , c is the local chord length of the blade element, in m, U ref is the local flow velocity of the blade element, in m / s; further, the actuating force per unit area on the rotor rotating actuating disk is calculated as: (9) In the formula, is the actuation force per unit area on the wind rotor rotating actuation disk, in N; B is the number of blades; is the area of ​​the annulus at different radial positions of the wind rotor, in m 2 ,in .

[0013] Furthermore, in step (4), the method for calculating the volume force of the wind turbine on the flow field is as follows: After obtaining the actuating force of the wind rotor surface, a three-dimensional Gaussian distribution function is used Smoothly distribute the force on the blade element to the surrounding grid and calculate the volume force of the wind turbine convection field , the calculation method is as follows: (10) (11) In the formula, is the volume force of the wind turbine convection field, in N; is the local grid volume, in m 3 ; For d is the three-dimensional Gaussian distribution function of the variable; d is the distance from the center of the grid to the blade element, in m; is the Gaussian smoothing parameter, which is used to adjust the concentration of the force, and its unit is m; According to Newton's third law, any cell in the area of ​​the actuator disk affected by the actuator disk p The driving force for: (12) Where N is the discrete number of the wind wheel in the radial direction; the negative sign indicates that the force on the fluid is opposite to the blade force.

[0014] Furthermore, in step (5), the velocity and density are solved as follows: (13) (14) In the formula, is the fluid density in kg / m 3 ; u is the macroscopic fluid velocity in m / s; m is the number of velocity concentration distribution functions in dimensionless units.

[0015] Here, the velocity of the fluid can be obtained by dividing the data obtained by formula (14) by the fluid density obtained by formula (13).

[0016] In the above, after the distribution function is obtained, the macroscopic physical quantity of the fluid is calculated by integrating the distribution function in the velocity space, and the macroscopic physical quantity at least includes the density and velocity of the fluid.

[0017] Secondly, the LBM equation can also be reduced to the macroscopic NS equation through the Chapman-Enskog method, so it can be used to solve macroscopic flow problems. The Chapman-Enskog method is the Chapman-Enskog method; Then, through the Chapman-Enskog method, the Lattice Boltzmann equation can be restored to the Navier-Stokes equation describing the macroscopic fluid motion, as shown in the following formula: (15) In the formula, is the viscosity coefficient, F is the external force term, ▽ p is the pressure gradient, and its relationship with the relaxation time in equation (5) is determined by the following formula: (16) Furthermore, step (1) of obtaining wind resource data is to obtain wind resource data of the wind farm area including at least wind speed and wind direction based on local meteorological data of the wind farm site or local mesoscale meteorological reanalysis data.

[0018] Furthermore, step (1) of obtaining wind resource data is to obtain wind resource data of the wind farm area including at least wind speed and wind direction based on local meteorological data of the wind farm site or local mesoscale meteorological reanalysis data.

[0019] Furthermore, in step (1), the method for constructing the wind farm computational domain is as follows: the location of the wind farm to be simulated is determined according to the longitude and latitude coordinates and a terrain geometry model is generated, on this basis an initial computational domain is formed, and the computational domain is discretized into a uniform Cartesian grid using an LBM unit type.

[0020] Furthermore, the LBM unit type is one of D3Q15, D3Q19 and D3Q27.

[0021] The LBM unit types include D3Q15, D3Q19, and D3Q27. Among them, D3Q15 refers to three-dimensional fifteen-speed, D3Q19 refers to three-dimensional nineteen-speed, and D3Q27 refers to three-dimensional twenty-seven-speed.

[0022] The present invention is beneficial in that: The present invention provides a wind farm wake simulation method based on LBM-LES. This method does not need to solve complex macroscopic flow motion equations, but instead solves the mesoscopic lattice Boltzmann equation, which greatly improves the solution efficiency compared to traditional methods. Secondly, the present invention adopts a rotating actuator disk method to simulate wind turbines, rather than using a superposition method to first calculate the flow field and then calculate the wind turbine wake, so that the wind farm wake simulation result is closer to the real wind field, and the calculation accuracy is significantly improved. BRIEF DESCRIPTION OF THE DRAWINGS

[0023] Figure 1 is a flow chart of the method of the present invention; Figure 2 is a schematic diagram of discrete blade elements; Figure 3 It is the velocity triangle and force diagram of the airfoil section; Figure 4 It is a schematic diagram of the rotating actuator disk model; Figure 5 The wind farm calculation domain and wind turbine layout diagram in Example 1; Figure 6 is the time-averaged wind speed profile curve at different positions from the wind rotor in Example 1; Figure 7 is a layout diagram of a wind farm in Example 2; Figure 8 is the time-averaged wind speed profile curve at different positions from the wind rotor in Example 2; Fig. 9 This is a diagram of wind farm power and thrust output in Example 2; Fig.10 The wind turbine power and thrust output time history in Example 2; Fig.11 This is a cloud diagram of the wind farm wake velocity at different times in Example 2. DETAILED DESCRIPTION

[0024] Example 1 This example is a comparison between the method described in the present invention and a specific experiment to verify the calculation accuracy of the method, as follows: This embodiment verifies the proposed LBM-LES wind farm wake efficient simulation method based on wind tunnel experimental measurement data. The wind tunnel experimental data used are the wind turbine model test results conducted in the low-speed backflow wind tunnel of the Norwegian University of Science and Technology at the 2011 Blind Test seminar. In the experiment, the size of the wind tunnel test section was 11.15 meters × 2.71 meters × 1.85 meters. The model used in the test is a three-blade wind turbine with a rotor diameter of 0.894 meters, installed on a tower with a height of 0.817 meters. The blades are manufactured with a NREL S826 airfoil with a relative thickness of 14% along the entire span, and the nacelle is a cylinder with a diameter of 90 mm, which is covered by a hemispherical hub in front and behind. The wind turbine is located at the centerline of the wind tunnel, and the rotor is 3.66 meters away from the entrance of the wind tunnel test section.

[0025] This study establishes a computational domain corresponding to the wind tunnel experiment, such as Figure 5 As shown: The inlet boundary condition is uniform flow, wind speed U ref =10m / s, air viscosity, the test wind turbine runs at the optimal tip speed ratio λ=6, at which the tip Reynolds number is, where is the tip chord length. The outlet boundary condition is zero gradient boundary, the surrounding boundaries are set as periodic conditions, and the time step is set to 2×10 -5 s.

[0026] Table 1 shows the comparison of the power coefficient and thrust coefficient of the BT1 wind turbine with the experimental data. It can be found from the table that the deviation of the power coefficient is 1.34%. Figure 6 The comparison between the normalized time-averaged wake velocity profile curves and the experimental measurement results at three positions downstream of the wind rotor plane - X=1D, X=3D, and X=5D - is shown. It can be found from the figure that in the near-wake area (X=1D), the calculated wake width is consistent with the wind tunnel experimental measurement results, which is about 1.1 times the diameter of the wind rotor. The above results show that the calculation accuracy of the proposed wind farm wake prediction method meets the requirements of engineering applications.

[0027] Table 1 Comparison of BT1 wind turbine power and thrust coefficient Lattice space Δ(m) <![CDATA[Grid number (×10 6 )]]> <![CDATA[ C P ]]> Medium grid 50 0.02 6.99 0.454 Wind tunnel test results - - - 0.448 Example 2 This embodiment applies the method of the present invention to a specific scenario to predict the wakes of multiple wind turbines; the details are as follows: Figure 7This is the layout diagram of a wind farm on flat terrain. The wind farm contains a total of 12 NREL 5MW wind turbines arranged in a 3×4 pattern. The wind turbines are arranged at a spacing of 5D in the incoming wind direction, where D=126m, which is the rotor diameter of the NREL 5MW. They are arranged at a spacing of 3D in the horizontal direction. The first row of wind turbines is 2D away from the inlet, and the last row is 8D away from the outlet. The computational domain is discretized using a D3Q19 uniform grid with Δx=2.5m, a length of 3.15km, or 25D, a width of 1.134km, or 9D, and a height of 0.378km, or 3D. The wind turbines operate at a rated wind speed of 11.4m / s and a rated speed of 12.1rpm. Based on the above setting parameters, the simulation calculation of the wake flow of the wind farm is carried out to obtain the aerodynamic force of the wind turbine rotor, as well as the macroscopic flow field density and velocity in the wind farm calculation domain. The Gaussian projection method is used to map the aerodynamic force to the volume force on the cell. The analysis is as follows: Figure 8 The figures in the figure are No. T2, No. T5, No. T8 and No. T11 wind turbines, respectively, showing the time-averaged wind speed profile curves at different positions of different rows of middle wind turbines, where four lines of different colors represent four wind turbines. It can be seen from the figure that the wake deficit of the first row of wind turbines, namely No. T2 wind turbine, is the most obvious, with the clearest wake structure, the longest propagation distance and the slowest wake recovery speed; however, for the second row of wind turbines, namely No. T5 wind turbine, due to the interference of the wake of the first row of wind turbines, the wake recovery is accelerated, resulting in a weakened wake deficit, a deformed wake structure and an unclear boundary; similarly, for the third row of wind turbines, namely No. T8 wind turbine, the wake interference effect of the first two rows of wind turbines further weakened the wake deficit; in the fourth row, namely No. T11 wind turbine, the comprehensive wake from the first three rows significantly increased the turbulence intensity of the inflow, resulting in a less obvious wake velocity profile, accompanied by a significant recovery of the average wake velocity. The wake disturbance also leads to the large-scale disruption of the wake vortex structure, forming small-scale vortices.

[0028] The power and thrust output of wind turbines in a wind farm are as follows: Fig. 9As shown in the figure, it can be clearly seen that the wake effect has a significant impact on the power and thrust output of the wind turbine. In terms of power output, the first row of wind turbines that are not affected by the wake reaches the rated power of 5.28MW; in contrast, the output of the second row of wind turbines is 1.21MW, which is only 22.49% of the output power of the first row of wind turbines. On the other hand, the power output of wind turbines located in the wake, such as the third and fourth rows of wind turbines, increased significantly, with the third row reaching about 2.1 MW and the fourth row recovering to about 3 MW. The wake interference between the upstream and downstream wind turbines accelerated the wake recovery behind the downstream wind turbine, resulting in a significant increase in the average wake velocity compared with the wake of a single wind turbine. Although the wake recovery was accelerated due to the wake interference, the wake effect still significantly affected the operation of the downstream wind turbines, and the power output of the fourth row of wind turbines still did not reach 60% of the power of the first row of wind turbines. Therefore, during the development and construction of the wind farm, a non-aligned layout, such as a staggered arrangement, should be adopted to minimize the wake impact on the downstream wind turbines. In the case where a linear layout is required due to terrain restrictions, increasing the distance between the first and second rows can provide sufficient time and space for the wake recovery of the first row of wind turbines, thereby increasing the total power generation of the wind farm.

[0029] Fig.10 The power-time curves of the wind turbines in the wind farm are shown. From the figure, we can see the effects of the speed deficit and wake interference caused by the wake on the power output. As the simulation time increases, the power output curve of the fourth row of wind turbines in the wind farm experiences three obvious bifurcations. When the wind farm starts to operate, all wind turbines gradually stabilize and their power converges to around the rated power of 5.28MW. However, around 50s of operation, the power output of the second, third and fourth rows of wind turbines drops sharply at the same time, dropping to about 1.2MW. This marks the first bifurcation of the power curve, because the strong wake deficit from the front row of wind turbines affects the downstream wind turbines, causing their power output to drop significantly. As the operating time continues to increase, around 125s, the power output curves of the third and fourth rows of wind turbines separate from the curves of the second row, which occurs when the interference wake generated by the wake interference reaches the third and fourth rows of wind turbines. In addition, at about 185 seconds of operation, the power curve experienced the last bifurcation, which was due to the mixed wake of the disturbed wind turbines reaching the fourth row of wind turbines, further accelerating the recovery of the wake velocity and increasing the average inflow velocity of the fourth row of wind turbines, thereby enhancing their power output. However, it is worth noting that although wake interference can accelerate wake recovery and enhance the power output of the wind farm, it also significantly increases the fluctuation of wind turbine power and structural vibration, which brings considerable pressure to the grid connection and daily maintenance of the wind farm.

[0030] Fig.11The transient wake velocity cloud diagram of the wind farm at different simulation times corresponds to the three key bifurcation points in the power-time curve: 50s, 125s and 185s; it can be found from the figure that in the initial stage of the wind farm operation, since all wind turbines experience uniform inflow, the wake structure is very similar. At 50s, the wake from the front row of wind turbines reaches the downstream wind turbines, and due to the wake effect, the inflow velocity of the downstream wind turbines is reduced. This causes the wind turbine power output to drop sharply, triggering the first bifurcation of the power output curve. As the operating time goes by, at 125s, the interference wakes formed by the first and second rows and the interference wakes formed by the second and third rows are very similar, reaching the third and fourth rows of wind turbines respectively, causing the average inflow velocity of the third and fourth rows of wind turbines to increase, triggering the second bifurcation of the power curve. When the operation time is further increased to 185s, the interference wakes from the first, second and third rows reach the fourth row, resulting in a change in the inflow conditions of the fourth row of wind turbines compared to the third row, which triggers the last bifurcation of the power curve. When the wind farm operation time reaches 300s, the inflow conditions of wind turbines in different rows have stabilized, and their velocity distribution is similar to that at 185s; the wake structure of the first row of wind turbines is the most obvious, the wake deficit of the second row of wind turbines is the largest, the wake deficit of the third row is partially recovered, and the recovery of the fourth row of wind turbines is obvious, and the wake deficit in the far wake area is the smallest.

Claims

1. A wind farm wake prediction method based on LBM-LES, characterized by: The following steps are included: (1) Obtain wind resources and wind farm terrain data and construct the wind farm calculation domain; (2) Applying boundary conditions and material parameters: Based on the data obtained from step (1), apply corresponding boundary conditions to the wind farm calculation domain and assign material parameters; (3) By replacing the effective viscosity in the collision model, the LES large eddy simulation model is embedded in the LBM, and the lattice Boltzmann equation is discretized and solved; (4) Calculate the aerodynamic force of the wind turbine rotor according to the wind turbine parameters in the wind farm, and use the Gaussian projection method to map the aerodynamic force into the volume force on the cell; (5) Calculate the macroscopic flow field density and velocity within the calculation domain of the wind farm and obtain the wind turbine wake velocity distribution information.

2. The method for predicting wind farm wake based on LBM-LES according to claim 1, characterized in that: The specific method of coupling the LES large eddy simulation model in LBM described in step (3) is as follows: The LES large eddy simulation model is coupled to the LBM through the sub-grid eddy viscosity model. The sub-grid eddy viscosity model is formed by superimposing the sub-grid viscosity coefficient considering the small-scale eddy dissipation on the original viscosity coefficient to form an equivalent turbulent viscosity. , thereby replacing the original viscosity coefficient and changing the collision model; （1） In the formula, is the original viscosity coefficient, in Pa·s, For The subgrid viscosity coefficient considering small-scale eddy dissipation is superimposed on the basis, and the unit is Pa·s; is the equivalent turbulent viscosity, in Pa·s; C s is the Smagorinsky constant, ranging from 0.1 to 0.2; is the filtering scale of the filter, taking the cube root value of the grid volume in LBM, that is, , Δ1, Δ2, Δ3 are the lengths of the unit cell in different directions, in m; defined by the following formula: （2） （3） In the formula, is the modulus of the strain rate tensor, is the strain rate tensor, in s -1 ; u i is the macroscopic fluid velocity, in m / s; i, j The value is 1, 2, or 3, corresponding to the x, y, and z directions respectively.

3. The method for predicting wind farm wake based on LBM-LES according to claim 1, characterized in that: Step (3) The discretization process of the lattice Boltzmann equation is: The continuous lattice Boltzmann equation, that is, Equation 4, is discretized to obtain the discrete form of the lattice Boltzmann equation, that is, Equation 5: （4） （5） In the formula, express t The velocity at position x is c i The number of microscopic particles is also called the distribution function; t is the time, the unit is s; c i is the velocity of microscopic particles, in m / s; F i is the force exerted by the wind turbine on the flow field, in N; ▽ is the gradient operator; As the collision operator, the Bhatnagar-Gross-Krook operator is used, that is, ; In the formula, is the local equilibrium distribution function in discrete velocity space is the weight coefficient defined in the selected discrete velocity model; is the fluid density in kg / m 3 ; u is the macroscopic fluid velocity, in m / s; c s is the lattice speed of sound, dimensionless; △ t is the time step, in seconds, is the relaxation time, dimensionless.

4. The method for predicting wind farm wake based on LBM-LES according to claim 3, characterized in that: In step (3), the lattice Boltzmann equation is evolved and solved using the relaxation migration method to obtain the distribution function after the update time step. The specific steps are as follows: First, the collision solution step is performed: （6） In the formula, The distribution function after solving for collisions; Then perform the migration evolution solution step: （7）。 5. The method for predicting wind farm wake based on LBM-LES according to claim 4, characterized in that: In step (4), the method for calculating the aerodynamic force of the wind turbine rotor is as follows: According to the momentum blade element theory, the blade is discretized into N blade element segments along the span direction. Assuming that the blade element segments do not interfere with each other, each blade element segment can be regarded as a two-dimensional airfoil extension. Furthermore, based on the incoming wind speed of the blade element segment and the aerodynamic parameters of the airfoil, the aerodynamic load of each blade element segment is calculated. The aerodynamic force of the blade element segment at any radial position r of the wind rotor is expressed as: （8） In the formula, F r is the aerodynamic force of the blade element at any radial position r of the wind rotor, in N / m; dr is the blade element length, in m; C L is the airfoil lift coefficient, C D is the airfoil drag coefficient, which can be obtained by looking up the table based on the airfoil information; e L is the direction of the lift unit vector, dimensionless, e D is the direction of the drag unit vector, and its unit is dimensionless; is the air density, in kg / m 3 , c is the local chord length of the blade element, in m, U ref is the local flow velocity of the blade element, in m / s; Furthermore, the actuation force per unit area on the wind rotor rotating actuation disk is calculated as: （9） In the formula, is the actuation force per unit area on the wind rotor rotating actuation disk, in N; B is the number of blades; is the area of ​​the annulus at different radial positions of the wind rotor, in m 2 ,in .

6. The method for predicting wind farm wake based on LBM-LES according to claim 5, characterized in that: In step (4), the method for calculating the volume force of the wind turbine on the flow field is as follows: After obtaining the actuating force of the wind rotor surface, a three-dimensional Gaussian distribution function is used Smoothly distribute the force on the blade element to the surrounding grid and calculate the volume force of the wind turbine convection field , the calculation method is as follows: （10） （11） In the formula, is the volume force of the wind turbine convection field, in N; is the local grid volume, in m 3 ; is a three-dimensional Gaussian distribution function with d as the variable; d is the distance from the center of the grid to the blade element, in m; is the Gaussian smoothing parameter, which is used to adjust the concentration of the force, and its unit is m; According to Newton's third law, any cell in the area of ​​the actuator disk affected by the actuator disk p The driving force for: （12） Where N is the discrete number of the wind wheel in the radial direction; the negative sign indicates that the force on the fluid is opposite to the blade force.

7. The method for predicting the wind farm wake based on LBM-LES according to claim 6, characterized in that: In step (5), the velocity and density are solved as follows: （13） （14） In the formula, is the fluid density in kg / m 3 ; u is the macroscopic fluid velocity in m / s; m is the number of velocity concentration distribution functions in dimensionless units.

8. The method for predicting wind farm wake based on LBM-LES according to claim 1, characterized in that: Step (1) of obtaining wind resource data is to obtain wind resource data of the wind farm area including at least wind speed and wind direction based on local meteorological data of the wind farm site or local mesoscale meteorological reanalysis data.

9. The method for predicting the wind farm wake based on LBM-LES according to claim 1, characterized in that: In step (1), the method for constructing the wind farm computational domain is as follows: the location of the wind farm to be simulated is determined according to the latitude and longitude coordinates and a terrain geometry model is generated, on this basis an initial computational domain is formed, and the wind farm computational domain is discretized into a uniform Cartesian grid using an LBM unit type.

10. The method for predicting the wind farm wake based on LBM-LES according to claim 9, characterized in that: The LBM unit type is one of D3Q15, D3Q19 and D3Q27.

Images