A wind farm wake modeling method considering atmospheric stratification and wind direction deflection

Abstract

The application discloses a wind farm wake modeling method considering atmospheric stratification and wind direction deflection, comprising the steps of obtaining basic parameters of a wind turbine and characteristic parameters of an atmospheric boundary layer; based on a coupling Ekman-surface atmospheric boundary layer analysis model, the atmospheric boundary layer is divided into an outer Ekman layer and an inner surface layer, and a flow velocity component, a spanwise velocity component, a boundary layer height, a friction velocity, a geostrophic wind component and a cross-isobar angle are calculated; based on a Townsend-Perry logarithmic scaling relationship, a flow turbulence intensity at a hub height is calculated, and an wake expansion rate is determined; a vortex sheet analysis method is used to construct a maximum speed loss amplitude of a yaw wake, a wake width and a yaw-induced displacement; a wind direction deflection displacement is calculated according to the flow velocity component and the spanwise velocity component, and is superimposed with the yaw-induced displacement to obtain a total wake displacement; and a spatial speed loss field is synthesized based on a Gaussian distribution to obtain a wind speed in a wake area and a power loss evaluation result.

Description

Technical Field

[0001] This invention relates to the field of wind farm wake modeling and atmospheric boundary layer analysis technology, and particularly to a wind farm wake modeling method that considers atmospheric stratification and wind direction deflection. Background Technology

[0002] When wind turbines are operating, they absorb kinetic energy from the incoming wind through their rotors, creating a wake region downstream where wind speed is reduced and turbulence is enhanced. When a downstream turbine is within the wake influence range of an upstream turbine, its effective inflow wind speed decreases, resulting in reduced turbine output power and consequently impacting the overall power generation efficiency and operational economy of the wind farm. Therefore, establishing a modeling method capable of rapidly calculating wake velocity distribution and wake propagation characteristics is crucial for wind farm layout design, turbine operation evaluation, and yaw control optimization.

[0003] Existing wind farm wake modeling methods mainly include analytical wake models and numerical simulation methods. Analytical wake models are fast and have strong engineering applicability, but they often rely on fixed empirical parameters and typically assume the wake to be an axisymmetric structure, making it difficult to reflect the influence of wind speed and direction variations with altitude in the atmospheric boundary layer on the wake propagation path. In particular, the Coriolis force causes wind direction to deflect with altitude within the boundary layer, resulting in lateral displacement and shear deformation of the wake during downstream propagation. Simultaneously, atmospheric thermal stratification alters turbulent mixing intensity and wake recovery rate, leading to differences in wake expansion patterns under stable, neutral, and unstable conditions. Existing analytical models, if they do not incorporate these factors into a unified calculation framework, are prone to errors in predicting the wake center location, velocity deficit magnitude, and downstream power loss.

[0004] Numerical simulation methods can describe complex flow field structures in detail, but they are computationally intensive and require high levels of mesh generation, boundary conditions, and computational resources, making them difficult to meet the needs of rapid assessment and online calculation for large-scale wind farms. With the increasing size of wind turbine generators, the increase in hub height, and the application of yaw control, the influence of Ekman spirals, stable stratification, and yaw-induced wake curl on wake distribution has become more prominent.

[0005] In view of this, the present invention provides a method for modeling wind farm wakes that takes into account atmospheric stratification and wind direction deflection. Summary of the Invention

[0006] The purpose of this invention is to provide a wind farm wake modeling method that considers atmospheric stratification and wind direction deflection, thereby solving the technical problem that existing analytical wake models cannot simultaneously couple atmospheric stratification, wind direction deflection, and yaw wake curling effects, resulting in insufficient accuracy in predicting wind farm wake velocity distribution and power loss.

[0007] To achieve the above objectives, the present invention provides the following technical solution: A method for modeling wind farm wakes that considers atmospheric stratification and wind direction deflection includes the following steps: Step S1: Obtain the basic parameters of the wind turbine and the characteristic parameters of the atmospheric boundary layer; Step S2: Based on the coupled Ekman-surface atmospheric boundary layer analysis model, the atmospheric boundary layer is divided into an outer Ekman layer and an inner surface layer. The velocity profile is coupled at the matching height to calculate the vertical distribution of the flow velocity component and the spanwise velocity component. The boundary layer height, friction velocity, geostrophic wind component and cross isobaric angle are output. Step S3: Calculate the turbulence intensity at the hub height based on the Townsend-Perry logarithmic scaling relationship, and determine the wake spread rate corresponding to atmospheric stability; Step S4: Calculate the maximum velocity deficit amplitude of the yaw wake, the polar angle-dependent wake width, and the spanwise displacement of the yaw-induced wake center based on the vortex plate analysis method. Step S5: Calculate the spanwise displacement of the wake of wind deflection based on the flow velocity component and the spanwise velocity component, and superimpose it with the spanwise displacement of the wake center induced by yaw to obtain the total wake displacement. Step S6: Based on the Gaussian distribution, the maximum velocity loss amplitude, the polar angle-related wake width, and the total wake displacement are synthesized into a spatial velocity loss field to obtain the wake region wind speed. The spatial velocity loss field is then applied to the downstream unit rotor sweep surface to obtain the equivalent inflow wind speed and power loss evaluation results. The wake region wind speed is obtained by subtracting the spatial velocity loss field from the background flow velocity that is not affected by the wake.

[0008] As a preferred technical solution of the present invention, the outer Ekman layer adopts the 3 / 2 power law total stress profile, and the inner surface layer adopts the Monin-Obukhov similarity theory. The coupled Ekman-surface atmospheric boundary layer analysis model introduces the Kazanskii-Monin stability parameter and the Zilitinkevich number to quantify the influence of surface cooling and free atmospheric thermal stratification on the atmospheric boundary layer structure, respectively. KazanskiiMonin stability parameters: ; Zilitinkevich number: ; in: The friction speed; It is von Kármán's constant; For Coriolis frequency, The BruntVäisälä frequency is for the free atmosphere.

[0009] As a preferred embodiment of the present invention, the boundary layer height is determined jointly by factors related to surface cooling, free atmospheric thermal stratification, and the neutral boundary layer, and satisfies the following: ; in: The dimensionless boundary layer height, ; The boundary layer height, For Coriolis frequency, For friction speed, For Kazanskii-Monin stability parameters, For Zilitinkevich numbers, , and The boundary layer height is an empirical constant. Substitute the Kazanskii-Monin stability parameters and Zilitinkevich numbers into the boundary layer height relationship and the geostrophic drag law, and iteratively solve for the geostrophic wind flow direction component, geostrophic wind spanwise component, friction velocity, and boundary layer height until the change in the corresponding parameters obtained from two adjacent iterations is less than the preset convergence threshold.

[0010] As a preferred technical solution of the present invention, the wake propagation rate is determined by the boundary layer height, hub height, friction speed, flow velocity at hub height and Townsend-Perry constant, so that the wake propagation rate changes with atmospheric stability and the wake propagation rate participates in the calculation of the polar angle-related wake width. Wherein, the flow velocity at the hub height is the value of the flow velocity component at the hub height, and the wake propagation rate adaptively changes with the boundary layer height, friction speed, or flow velocity at the hub height.

[0011] As a preferred embodiment of the present invention, the polar angle-dependent wake width is determined by the wake propagation rate, the downstream distance, and the wake shape function. The wake shape function is obtained by multiplying the initial shape function and the dimensionless vortex plate position. The downstream distance is the distance between the location to be calculated and the center of the upstream wind turbine rotor along the incoming flow direction. The polar angle-dependent wake width is used to characterize the local width of the wake boundary at the same downstream location as the polar angle changes.

[0012] As a preferred technical solution of the present invention, the initial shape function is determined by the wind turbine radius, the wind turbine effective wind-receiving area correction coefficient, the yaw angle and the polar angle; the position of the dimensionless vortex plate is calculated by the analytical expression of the truncated power series expansion to characterize the non-uniform expansion shape of the yaw wake in different polar angle directions. The wind turbine's effective wind-receiving area correction coefficient is used to correct the change in the equivalent wind-receiving area caused by the yaw angle, and the position of the dimensionless vortex plate changes with the dimensionless wake evolution time and polar angle. The dimensionless wake evolution time is determined by the wind turbine's effective wind-receiving area correction coefficient, thrust coefficient, yaw angle, friction speed, flow velocity at hub height, flow velocity at the height to be calculated, downstream distance, and wind turbine radius.

[0013] As a preferred embodiment of the present invention, the maximum velocity loss amplitude is determined according to the segmentation of the wake development region; When the downstream calculation location is in the potential core region, the maximum velocity loss amplitude is determined according to the velocity loss amplitude in the potential core region corresponding to the axial induction factor. When the downstream calculation location is in the attenuation zone, the maximum velocity loss amplitude is determined based on the thrust coefficient, yaw angle, polar angle-related wake width at the corresponding downstream location, and rotor radius. The transition position between the potential core region and the decay region is determined based on the continuous matching relationship between the velocity loss amplitude of the potential core region and the velocity loss amplitude of the decay region. The transition position is used to ensure that the maximum velocity loss amplitude transitions continuously between the potential core region and the decay region.

[0014] As a preferred embodiment of the present invention, the spanwise displacement of the wake caused by wind deflection... Satisfy the following formula: ; in, downstream distance, The height of the position to be calculated. For height The flow velocity component at that location, For height The spanwise velocity component at that location; The total wake displacement satisfies: ; in, This represents the total wake displacement. This refers to the spanwise displacement of the wake caused by wind deflection. For the yaw-induced wake center spanwise displacement, the wake spanwise displacement caused by wind deflection and the yaw-induced wake center spanwise displacement use the same positive spanwise coordinate direction.

[0015] As a preferred embodiment of the present invention, the wind speed at any location within the wake region is obtained by subtracting the spatial velocity deficit field from the streamwise velocity component within the atmospheric boundary layer, and satisfies the following: ; in, Location within the wake region Wind speed at the location, For spatial velocity deficit fields; The spatial velocity deficit field determines the spanwise position of the wake center using the total wake displacement, determines the range of wake expansion in different polar angle directions using the polar angle-related wake width, and determines the degree of velocity deficit in the wake center region using the maximum velocity deficit amplitude.

[0016] As a preferred technical solution of the present invention, after completing the calculation of boundary layer velocity profile, wake spread rate, yaw wake parameters, wind direction deflection displacement, total wake displacement and spatial velocity deficit field, the wind speed, spatial velocity deficit field, maximum velocity deficit amplitude, wake width, total wake displacement and power loss evaluation results at any location in the wake region are output. The downstream wind turbine rotor swept surface is divided into multiple calculation units. The corresponding local flow wind speed is determined based on the position of each calculation unit in the spatial velocity deficit field. The local flow wind speed of each calculation unit is then weighted by area to obtain the equivalent inflow wind speed of the downstream wind turbine. The equivalent inflow wind speed is input into the wind turbine power curve to obtain the output power under the wake influence condition. The power loss evaluation result is determined based on the difference between the output power and the reference output power when not affected by the wake.

[0017] Compared with the prior art, the present invention has the following beneficial effects: This invention obtains a boundary layer velocity profile containing both streamwise and spanwise velocity components by coupling the Ekman-surface atmospheric boundary layer analysis model, enabling the wake modeling process to incorporate the effects of atmospheric stratification and wind deflection caused by the Coriolis force. It establishes a correspondence between wake propagation rate and atmospheric stability using the Townsend-Perry logarithmic scaling relationship, freeing the wake recovery process from reliance on fixed empirical parameters. Furthermore, it calculates the maximum velocity deficit amplitude, wake width, and yaw-induced displacement using a vortex-plate-based yaw wake model, and synthesizes the spatial velocity deficit field by superimposing the wind deflection displacement and the yaw-induced displacement. This allows for the acquisition of wake velocity distribution and power loss evaluation results under different yaw angles and atmospheric boundary layer conditions, improving the applicability of wind farm wake prediction, layout optimization, and operation assessment. Attached Figure Description

[0018] To more clearly illustrate the technical solutions in the embodiments of this application or the prior art, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments recorded in this invention. For those skilled in the art, other drawings can be obtained based on these drawings.

[0019] Figure 1 This is a flowchart of a wind farm wake modeling method that considers atmospheric stratification and wind direction deflection according to the present invention. Figure 2 This is a comparison of wake velocity deficit distribution at different downstream distances x / D under neutral boundary layer conditions and yaw angle β = -20° in an embodiment of the present invention. Figure 3 This is a comparison diagram of wake velocity deficit distribution at different downstream distances x / D under neutral boundary layer conditions and yaw angle β=0° in an embodiment of the present invention. Figure 4 This is a comparison diagram of wake velocity deficit distribution at different downstream distances x / D under neutral boundary layer conditions and yaw angle β=20° in an embodiment of the present invention. Detailed Implementation

[0020] The technical solution of the present invention will be described below with reference to the accompanying drawings. It should be understood that the following embodiments are used to illustrate the implementation of the present invention and are not intended to limit the scope of protection of the present invention. Equivalent substitutions made by those skilled in the art regarding parameter values, solution accuracy, mesh scale, or multi-machine overlay methods without departing from the technical concept of the present invention should all fall within the scope of protection of the present invention.

[0021] Example 1 like Figure 1 As shown, this embodiment provides a wind farm wake modeling method that considers atmospheric stratification and wind direction deflection. By coupling the Ekman-surface atmospheric boundary layer analysis model with a vortex-plate-based yaw wake model, the atmospheric boundary layer velocity profile, wind direction deflection correction, wake propagation rate, yaw-induced wake curling, and Gaussian velocity deficit field are correlated and calculated to obtain the downstream wake velocity deficit distribution and power loss evaluation results of the wind turbine. The method includes the following steps: Step S101: Obtain the basic parameters of the wind turbine and the characteristic parameters of the atmospheric boundary layer.

[0022] The basic parameters of the wind turbine include rotor radius, rotor diameter, hub height, thrust coefficient, axial induction factor, and yaw angle. These basic parameters are used to calculate the maximum velocity loss amplitude of the yaw wake, the polar angle-dependent wake width, and the spanwise displacement of the yaw-induced wake center.

[0023] The atmospheric boundary layer characteristic parameters include surface roughness, Coriolis frequency, Kazanskii-Monin stability parameter, and Zilitinkevich number. These atmospheric boundary layer characteristic parameters are used to calculate the boundary layer height, friction velocity, geostrophic wind component, and the vertical distribution of the flow velocity component and spanwise velocity component.

[0024] In specific calculations, the incoming flow direction is used as the flow direction, the horizontal direction perpendicular to the incoming flow direction is used as the spanwise direction, and the vertical height direction is used as the height direction. A three-dimensional coordinate system is established with the wind turbine center as the coordinate reference, where x represents the downstream flow coordinate, y represents the spanwise coordinate, z represents the vertical height above the ground, and the position of the wind turbine center is represented as... Under this coordinate relationship, the position in the downstream wake region of the wind turbine is calculated.

[0025] Step S102: Calculate the atmospheric boundary layer velocity profile based on the coupled Ekman-surface atmospheric boundary layer analysis model.

[0026] Specifically, the atmospheric boundary layer is divided into the outer Ekman layer and the inner surface layer. The outer Ekman layer is described by a 3 / 2 power-law total stress profile, which describes the attenuation of turbulent stress with height; that is, the total stress amplitude satisfies: ; in: The total surface stress is given by , and h is the boundary layer height. The near-surface velocity profile of the inner surface layer is described using the Monin-Obukhov similarity theory, and its flow velocity gradient satisfies: ; in: For friction speed, It is von Kármán's constant. The length of Obukhov, This is a universal stability function. The two layers are matched at height z. m Velocity profile coupling is performed at 0.2h to ensure that the velocities and velocity gradients of the two layers are continuous at the matching height, thereby avoiding abrupt changes in velocity values ​​or velocity gradients between the outer Ekman layer and the inner surface layer at the matching position.

[0027] The model incorporates the KazanskiiMonin stability parameter. and Zilitinkevich numbers The effects of surface cooling and free atmospheric thermal stratification on the atmospheric boundary layer structure were quantified separately.

[0028] ; ; in: The friction speed; It is von Kármán's constant; For Coriolis frequency, BruntVäisälä frequency for free atmosphere; according to and The model can adaptively describe the differences in velocity profiles under stable, neutral, and unstable boundary layer conditions by taking the value of .

[0029] Boundary layer height h Determined jointly by factors related to surface cooling, free atmospheric thermal stratification, and the neutral boundary layer, the following relationship is satisfied: ; in: The dimensionless boundary layer height, ; The boundary layer height, For Coriolis frequency, For friction speed, For KazanskiiMonin stability parameters, For Zilitinkevich numbers, , and This is an empirical constant for the boundary layer height.

[0030] The geostrophic wind component and friction velocity are solved using the Geostrophic Drag Law (GDL). Given the geostrophic wind amplitude G, the geostrophic wind flow direction component... U g Geospin wind direction component V g and friction speed u * Solve the following GDL equations iteratively: ; in z 0 For surface roughness, A and B The coefficients, which depend on stability parameters, boundary layer height, and friction velocity, are specifically expressed by the coupling conditions between the outer Ekman layer and the inner surface layer at the matching height. During the iteration process, the following is first given... u * and h Initial value, calculate A , B and ĥ Then update via GDL equations u * Then update through boundary layer height relationship h Iterate repeatedly until u * , U g , V g and hBoth convergences are defined as the change in the corresponding parameters obtained in two adjacent iterations being less than a preset convergence threshold.

[0031] After iterative convergence, the obtained u * , h , U g , V g Substituting the velocity profile expressions for the outer Ekman layer and the inner surface layer, we can obtain the complete streamwise velocity components within the atmospheric boundary layer. U(z) and spanwise velocity components V(z) The vertical distribution of the boundary layer is also output. h Friction speed u * Geostrophic wind direction component U g Geospin wind direction component V g and cross isobars α 0 =tan -1 (V g / U g ) .

[0032] Through the above processing, the output results form two subsequent calculation paths: one path will... U(z) , u * , h and U hub Used in step S103 to calculate the turbulence intensity and wake propagation rate; another path will U(z) and V(z) Used in step S105 to calculate the spanwise displacement of the wake caused by wind deflection.

[0033] Step S103: Calculate the directional turbulence intensity at the hub height based on the Townsend-Perry logarithmic scaling relation, and establish the functional relationship between the wake spread rate kW and atmospheric stability.

[0034] In this step, the flow-direction turbulence intensity at the hub height is calculated based on the Townsend-Perry logarithmic scaling relation, and the wake spread rate is established based on this flow-direction turbulence intensity. The adaptive function relationship between atmospheric stability and atmospheric stability: ; in, The Townsend-Perry constant is set to 1.25. , n=6; This refers to the hub height of the wind turbine. Let be the flow velocity at the height of the wheel hub.

[0035] Through the above processing, the wake spread rate is... The wake propagation rate is determined by the atmospheric boundary layer structure and the turbulence characteristics at the hub height, thus relating the wake propagation process to atmospheric stability. When the atmospheric stability changes and causes changes in the boundary layer height, friction velocity, or flow velocity at the hub height, the wake propagation rate kwk_wkw changes accordingly, rather than using a fixed empirical value under different operating conditions.

[0036] Step S104: Using a vortex plate-based analysis method, an analytical model of the yaw wind turbine wake is constructed to calculate the maximum velocity deficit magnitude C(x) and the polar angle-dependent wake width of the yaw wind turbine wake. and the spanwise displacement of the wake center ŷ c The parsing expression for (ŧ).

[0037] The wake width is expressed as: ; in: downstream distance, Polar angle, This is the wake shape function. The polar angle-dependent wake width refers to the wake width at the same flow direction position. The wake boundary at the polar angle varies with the polar angle The varying local width, hereinafter referred to as wake width, refers to the wake width related to the polar angle.

[0038] The wake shape function is determined by the initial shape function and the position of the dimensionless vortex plate; Wake shape function: ; Initial shape function: : Among them: the position of the dimensionless vortex plate The non-uniform expansion morphology of the yaw wake is characterized by calculating the analytical expression of the truncated power series expansion in different polar angle directions. Through this process, the wake width is no longer restricted to an axisymmetric distribution, but can reflect the wake curling and lateral non-uniform expansion characteristics under yaw operation.

[0039] The spanwise displacement of the wake center induced by yaw is determined by the lateral velocity induced by the vortex plates. Define the dimensionless wake evolution time. in This is the correction factor for the effective wind-receiving area of ​​the wind turbine. For thrust coefficient, Let be the rotor radius. This dimensionless time characterizes the evolution of the vortex plate under the influence of turbulent dissipation and mean flow shear during downstream propagation. Based on this dimensionless time, the spanwise displacement of the wake center induced by yaw is calculated. The parsing expression is: The first term is the mainstream lateral induced displacement, and the second term is the correction caused by the ground-mirror vortex. Therefore, the spanwise displacement of the yaw-induced wake center is denoted as... Its dimensions are restored to: ; in, The sign of the yaw is determined according to the positive direction of the spanwise coordinate established in step S101, so that the displacement direction remains consistent under positive and negative yaw conditions.

[0040] The maximum velocity deficit magnitude C(x) is determined segmentally according to the wake development region. When the downstream calculation location is in the potential core region, the maximum velocity deficit magnitude is determined according to the velocity deficit magnitude in the potential core region corresponding to the axial induction factor. In the core region hour, ; In the attenuation region hour, ; in, This represents the transition point between the potential core region and the decay region. This represents the wake width at the corresponding downstream location. The transition point between the potential core region and the attenuation region. The maximum velocity loss amplitude is determined based on the continuous matching relationship between the velocity loss amplitude in the potential core region and the velocity loss amplitude in the decay region, so that the maximum velocity loss amplitude transitions continuously between the potential core region and the decay region.

[0041] The above calculations yielded the velocity loss magnitude, wake width, and yaw-induced spanwise displacement of the yaw wake at different downstream locations.

[0042] Step S105: Calculate the spanwise displacement of the wake caused by wind deflection and obtain the total wake displacement. Based on the flow velocity component U(z) and spanwise velocity component V(z) obtained in step S102, calculate the spanwise displacement of the wake caused by wind deflection: ; in, This refers to the spanwise displacement of the wake caused by wind deflection. downstream distance, The height of the position to be calculated. For height The flow velocity component at that location, For height The spanwise velocity component at that point; this calculation is equivalent to calculating based on the height. The ratio of the spanwise velocity to the flow velocity at a given point determines the direction of wake propagation, and the downstream propagation distance is also considered. Determine the cumulative spanwise displacement caused by wind deflection.

[0043] Subsequently, the wake caused by the wind deflection will be shifted longitudinally. With yaw-induced displacement Linear superposition yields the total wake displacement. ; ; in, , and All displacements are calculated using the same positive spanwise coordinate determined in step S101. When two displacements are in the same direction, their absolute values ​​are added together; when two displacements are in opposite directions, they are algebraically superimposed based on their signs. Wind deflection correction, as an independent branch directly derived from the boundary layer velocity profile, merges with yaw-induced displacement in the total wake displacement calculation stage, enabling the wind deflection caused by the Coriolis force to enter the subsequent velocity deficit field synthesis process.

[0044] Step S106: Based on the Gaussian distribution assumption, synthesize the spatial velocity loss field and calculate the wind speed in the wake region. Based on the Gaussian distribution assumption, the maximum velocity loss amplitude C(x) and the wake width are... and total wake displacement Synthesized into a spatial velocity deficit field This allows us to obtain the wind speed at any location within the wake region: ; in: Location within the wake region Wind speed at the location, The background flow velocity at height zzz, unaffected by the wake. This is a space velocity deficit field.

[0045] In the process of synthesizing the spatial velocity deficit field, the total wake displacement is used as the basis. Determine the position of the wake center in the spanwise direction at the corresponding height, and correlate the wake width with the polar angle. The extent of the wake's extension in different polar angle directions is determined, and the degree of velocity deficit in the wake's central region is determined by the maximum velocity deficit amplitude C(x), thereby forming a three-dimensional velocity deficit field that varies with the flow direction, spanwise position, and height.

[0046] To evaluate the power loss of downstream wind turbines, the swept surface of the downstream wind turbine rotor is divided into multiple computational units. Based on the position coordinates of the center of each computational unit, the corresponding local flow velocity is obtained from the spatial velocity deficit field. Then, the local flow velocity of each computational unit within the swept surface is area-weighted to obtain the equivalent inflow velocity of the downstream wind turbine. This equivalent inflow velocity is compared with the reference inflow velocity when not affected by the wake of the upstream turbine, and combined with the power curve of the wind turbine, the output power under wake influence conditions and the power loss evaluation results relative to the reference condition are obtained.

[0047] When the method is used in a wind farm consisting of multiple wind turbine units, the spatial velocity loss field of each upstream wind turbine unit on the rotor sweep surface of the target downstream wind turbine unit is determined sequentially according to the direction of the incoming flow, and the velocity loss influence of multiple upstream wakes is synthesized by using a preset wake superposition method.

[0048] Therefore, the calculations of boundary layer velocity profile, wake expansion rate, yaw wake parameters, wind direction deflection displacement, total wake displacement, and spatial velocity deficit field are completed. The results output include wind speed, spatial velocity deficit field, maximum velocity deficit magnitude, wake width, total wake displacement, and power loss evaluation at any location within the wake region. By treating the atmospheric boundary layer velocity profile, wind direction deflection displacement, yaw-induced displacement, and velocity deficit field as a continuous calculation process, the modeling and evaluation of wind farm wake distribution and power loss are achieved.

[0049] Example 2 This embodiment, based on Embodiment 1, illustrates the wake velocity deficit distribution under different yaw angle conditions in the neutral boundary layer (CNBL) environment.

[0050] The simulation conditions used in this embodiment are: simulating conventional neutral boundary layer (CNBL) atmospheric conditions, with a simulation domain size of: Grid resolution: Fan rotor diameter: Wheel hub height: Earthslide wind amplitude: Coriolis frequency Surface roughness The parameters for each operating condition are summarized in Table 1.

[0051] Table 1: Operating Parameters Under the above conditions, yaw angles β = -20°, β = 0°, and β = 20° were set respectively, and the normalized wake velocity deficit distribution was analyzed at different downstream distances x / D. A comparison will be made; specifically: ; in: For height The flow velocity unaffected by the wake. Location within the wake region The flow velocity at that location, The velocity of the incoming flow at the height of the wheel hub.

[0052] Figures 2 to 4 In the diagram, different columns correspond to different downstream distances x / D, and different rows correspond to different correction methods. "Yaw and Wind Direction Deflection" indicates that both yaw-induced correction and wind direction deflection correction are used simultaneously; "Yaw" indicates that only yaw-induced correction is used without wind direction deflection correction; and "No Correction" indicates that neither yaw-induced correction nor wind direction deflection correction is used. The horizontal axis in the diagram is the spanwise coordinate y / D, and the vertical axis is the relative hub height. .

[0053] The operating condition with a yaw angle β = -20°, such as Figure 2 As shown, under neutral boundary layer conditions and under "uncorrected" conditions, the velocity loss contour lines at each downstream distance are centered at the hub (y=0, The wake is centrally symmetrically distributed, with the loss amplitude decreasing as x / D increases and the width gradually increasing. Under "yaw" conditions, the wake center shifts in the positive y direction (due to the positive lateral force generated by negative yaw), and the contour lines exhibit a slight "C"-shaped curl. This shift direction is determined by the yaw angle notation convention used in this embodiment and the positive direction of the spanwise coordinates.

[0054] Under the conditions of "yaw and wind deflection", in addition to the positive y offset caused by yaw, wind deflection (V(z) in the boundary layer changes with altitude) causes different lateral displacements in the wake at different altitudes, resulting in the contour lines being twisted in the vertical direction, exhibiting shear deformation, and the maximum loss point is no longer located on the same vertical line, and the wake as a whole is tilted in the positive y direction.

[0055] It can be seen that under negative yaw conditions, yaw-induced displacement mainly changes the overall spanwise position of the wake center, while wind deflection displacement further introduces differential displacement that varies with altitude; the superposition of the two changes both the overall position and vertical shear pattern of the wake.

[0056] Under the condition of yaw angle β = 0°, such as Figure 3As shown, under neutral boundary layer conditions, when the yaw angle β = 0°, the results of the "yaw" condition and the "no correction" condition almost coincide (because the yaw-induced force is zero), and the wake maintains an axisymmetric shape. The "yaw and wind deflection" condition, however, separately shows the effect of wind deflection: since V(z) is close to zero near the hub height, but negative above the hub (in the Northern Hemisphere) and positive below, the upper half of the wake deflects in the negative y direction, and the lower half deflects in the positive y direction, exhibiting an overall "S"-shaped shear deformation, and the amount of deflection accumulates as x / D increases.

[0057] Therefore, the zero-yaw condition can be used to separate and observe the independent effect of wind deflection correction on wake distribution, avoiding interference from yaw-induced displacement on the analysis results.

[0058] Under the condition of yaw angle β = 20°, such as Figure 4 As shown, under neutral boundary layer conditions, positive yaw generates a lateral force in the negative y direction. Therefore, under yaw conditions, the wake center shifts in the negative y direction, and the curling direction is the same as... Figure 2 Conversely, under the conditions of "yaw and wind deflection," the directions of wind deflection and yaw-induced displacement are the same in the upper half of the hub (negative y-axis), and opposite in the lower half (yaw towards negative y-axis, wind deflection towards positive y-axis). Therefore, the upper half of the wake is significantly offset greater than the lower half, and the contour lines are strongly distorted, forming a distinct asymmetric "hook-like" structure. The comparison of the three lines shows that neither yaw nor wind deflection alone can accurately describe this complex deformation; only by considering both simultaneously can the asymmetric curling and shear coupling morphology observed in LES or actual measurements be reproduced.

[0059] This embodiment demonstrates that, through Figures 2 to 4 In comparison, under neutral boundary layer conditions, the wake velocity deficit distribution differs for different yaw angle conditions. By combining the Ekman surface atmospheric boundary layer analysis model, the TownsendPerry logarithmic scaling relation, the vortex plate yaw wake model, and the wind deflection correction term, the effects of atmospheric stratification, wind deflection, and yaw wake curling on the wake velocity deficit field can be described in the same modeling process.

[0060] This embodiment demonstrates that, through Figures 2 to 4 In comparison, under neutral boundary layer conditions, the wake velocity deficit distribution differs for different yaw angle conditions. By combining the Ekman-surface atmospheric boundary layer analysis model, the Townsend-Perry logarithmic scaling relation, the vortex plate yaw wake model, and the wind deflection correction term, the effects of atmospheric stratification, wind deflection, and yaw wake curling on the wake velocity deficit field can be described in the same modeling process.

[0061] The foregoing has shown and described the basic method, main features, and advantages of the present invention. Those skilled in the art should understand that the present invention is not limited to the above examples; the examples and descriptions in the specification are merely illustrative of the method. Various changes and modifications can be made to the present invention without departing from its spirit and scope, and all such changes and modifications fall within the scope of the present invention as claimed. The scope of protection of this invention is defined by the appended claims and their equivalents.

Claims

1. A method for modeling wind farm wakes that considers atmospheric stratification and wind direction deflection, characterized in that, Includes the following steps: Step S1: Obtain the basic parameters of the wind turbine and the characteristic parameters of the atmospheric boundary layer; Step S2: Based on the coupled Ekman-surface atmospheric boundary layer analysis model, the atmospheric boundary layer is divided into an outer Ekman layer and an inner surface layer. The velocity profile is coupled at the matching height to calculate the vertical distribution of the flow velocity component and the spanwise velocity component. The boundary layer height, friction velocity, geostrophic wind component and cross isobaric angle are output. Step S3: Calculate the turbulence intensity at the hub height based on the Townsend-Perry logarithmic scaling relationship, and determine the wake spread rate corresponding to the atmospheric stability. Step S4: Calculate the maximum velocity deficit amplitude of the yaw wake, the polar angle-dependent wake width, and the spanwise displacement of the yaw-induced wake center based on the vortex plate analysis method. Step S5: Calculate the spanwise displacement of the wind direction deflection wake based on the flow velocity component and the spanwise velocity component, and superimpose it with the spanwise displacement of the yaw-induced wake center according to the unified positive direction of the spanwise coordinate to obtain the total wake displacement. Step S6: Based on the Gaussian distribution, the maximum velocity loss amplitude, the polar angle-related wake width, and the total wake displacement are synthesized into a spatial velocity loss field to obtain the wake region wind speed. The spatial velocity loss field is then applied to the downstream unit rotor sweep surface to obtain the equivalent inflow wind speed and power loss evaluation results. The wake region wind speed is obtained by subtracting the spatial velocity loss field from the background flow velocity that is not affected by the wake.

2. The wind farm wake modeling method considering atmospheric stratification and wind direction deflection according to claim 1, characterized in that, The outer Ekman layer adopts a 3 / 2 power-law total stress profile, and the inner surface layer adopts the Monin-Obukhov similarity theory. The coupled Ekman-surface atmospheric boundary layer analysis model introduces the Kazanskii-Monin stability parameter and the Zilitinkevich number to quantify the influence of surface cooling and free atmospheric thermal stratification on the atmospheric boundary layer structure, and determines the velocity profile differences under stable, neutral or unstable boundary layer conditions based on the Kazanskii-Monin stability parameter and the Zilitinkevich number. KazanskiiMonin stability parameters: ; Zilitinkevich number: ; in: The friction speed; It is von Kármán's constant; For Coriolis frequency, The BruntVäisälä frequency is for the free atmosphere.

3. The wind farm wake modeling method considering atmospheric stratification and wind direction deflection according to claim 2, characterized in that, The boundary layer height is determined by a combination of factors related to surface cooling, free atmospheric thermal stratification, and the neutral boundary layer, and satisfies the following: ； in: The dimensionless boundary layer height. ; The boundary layer height, For Coriolis frequency, For friction speed, For Kazanskii-Monin stability parameters, For Zilitinkevich numbers, , and The boundary layer height is an empirical constant. Substitute the Kazanskii-Monin stability parameters and Zilitinkevich numbers into the boundary layer height relationship and the geostrophic drag law, and iteratively solve for the geostrophic wind flow direction component, geostrophic wind spanwise component, friction velocity, and boundary layer height until the change in the corresponding parameters obtained from two adjacent iterations is less than the preset convergence threshold.

4. The wind farm wake modeling method considering atmospheric stratification and wind direction deflection according to claim 1, characterized in that, The wake propagation rate is determined by the boundary layer height, hub height, friction speed, flow velocity at hub height, and Townsend-Perry constant, so that the wake propagation rate changes with atmospheric stability and participates in the calculation of the polar angle-dependent wake width. Wherein, the flow velocity at the hub height is the value of the flow velocity component at the hub height, and the wake propagation rate adaptively changes with the boundary layer height, friction speed, or flow velocity at the hub height.

5. The wind farm wake modeling method considering atmospheric stratification and wind direction deflection according to claim 1, characterized in that, The polar angle-dependent wake width is determined by the wake propagation rate, downstream distance, and wake shape function. The wake shape function is obtained by multiplying the initial shape function and the dimensionless vortex plate position. The downstream distance is the distance between the location to be calculated and the center of the upstream wind turbine rotor along the incoming flow direction. The polar angle-dependent wake width is used to characterize the local width of the wake boundary at the same downstream location as the polar angle changes.

6. The wind farm wake modeling method considering atmospheric stratification and wind direction deflection according to claim 5, characterized in that, The initial shape function is determined by the rotor radius, the effective wind-receiving area correction coefficient of the rotor, the yaw angle, and the polar angle; the position of the dimensionless vortex plate is calculated by the analytical expression of the truncated power series expansion to characterize the non-uniform expansion shape of the yaw wake in different polar angle directions. The wind turbine's effective wind-receiving area correction coefficient is used to correct the change in the equivalent wind-receiving area caused by the yaw angle, and the position of the dimensionless vortex plate changes with the dimensionless wake evolution time and polar angle. The dimensionless wake evolution time is determined by the wind turbine's effective wind-receiving area correction coefficient, thrust coefficient, yaw angle, friction speed, flow velocity at hub height, flow velocity at the height to be calculated, downstream distance, and wind turbine radius.

7. The wind farm wake modeling method considering atmospheric stratification and wind direction deflection according to claim 1, characterized in that, The maximum velocity loss amplitude is determined according to the segmentation of the wake development region; When the downstream calculation location is in the potential core region, the maximum velocity loss amplitude is determined according to the velocity loss amplitude in the potential core region corresponding to the axial induction factor. When the downstream calculation location is in the attenuation zone, the maximum velocity loss amplitude is determined based on the thrust coefficient, yaw angle, polar angle-related wake width at the corresponding downstream location, and rotor radius. The transition position between the potential core region and the decay region is determined based on the continuous matching relationship between the velocity loss amplitude of the potential core region and the velocity loss amplitude of the decay region. The transition position is used to ensure that the maximum velocity loss amplitude transitions continuously between the potential core region and the decay region.

8. The wind farm wake modeling method considering atmospheric stratification and wind direction deflection according to claim 1, characterized in that, Spanwise displacement of wake caused by wind deflection Satisfy the following formula: ; in, downstream distance, The height of the position to be calculated. For height The flow velocity component at that location, For height The spanwise velocity component at that location; The total wake displacement satisfies: ; in, This represents the total wake displacement. This refers to the spanwise displacement of the wake caused by wind deflection. For the yaw-induced wake center spanwise displacement, the wake spanwise displacement caused by wind deflection and the yaw-induced wake center spanwise displacement use the same positive spanwise coordinate direction.

9. The wind farm wake modeling method considering atmospheric stratification and wind direction deflection according to claim 1, characterized in that, The wind speed at any location within the wake region is obtained by subtracting the spatial velocity deficit field from the streamwise velocity component within the atmospheric boundary layer, and satisfies: ； in, Location within the wake region Wind speed at the location, This is a space velocity deficit field; The spatial velocity deficit field determines the spanwise position of the wake center using the total wake displacement, determines the range of wake expansion in different polar angle directions using the polar angle-related wake width, and determines the degree of velocity deficit in the wake center region using the maximum velocity deficit amplitude.

10. A wind farm wake modeling method considering atmospheric stratification and wind direction deflection according to claim 1, characterized in that, After calculating the boundary layer velocity profile, wake spread rate, yaw wake parameters, wind deflection displacement, total wake displacement, and spatial velocity deficit field, the system outputs the wind speed, spatial velocity deficit field, maximum velocity deficit magnitude, wake width, total wake displacement, and power loss evaluation results at any location within the wake region. The downstream wind turbine rotor swept surface is divided into multiple calculation units. The corresponding local flow wind speed is determined based on the position of each calculation unit in the spatial velocity deficit field. The local flow wind speed of each calculation unit is then weighted by area to obtain the equivalent inflow wind speed of the downstream wind turbine. The equivalent inflow wind speed is input into the wind turbine power curve to obtain the output power under the wake influence condition. The power loss evaluation result is determined based on the difference between the output power and the reference output power when not affected by the wake.

Images