Novel yaw fan wake flow calculation model method

Abstract

The invention provides a novel yawing fan wake flow calculation model method, and relates to the technical field of fan speed field calculation, a near wake flow field of a fan is simulated and calculated by using a CFD method based on an actuating line theory, then a wake flow center is determined by using a zero point searching method, and then a far wake flow field of the fan is simulated and calculated by using an improved Jensen model. And smoothing calculation is carried out on the coupling region by adopting an exponential factor to neutralize CFD and Jensen. Compared with a medium-precision wake flow model, the coupling method is higher in calculation precision; compared with a high-precision model, the coupling method is lower in calculation cost and is more suitable for preliminary estimation calculation in engineering.

Description

Technical Field

[0001] The present invention relates to the technical field of fan speed field calculation, and in particular to a new yaw fan wake calculation model method. Background Art

[0002] With the continuous progress of wind power generation technology, the scale of wind farms is becoming increasingly large, the mutual influence between wind turbines is becoming increasingly complex, and the mutual influence between wind farms is also growing. In order to pursue higher energy conversion efficiency, modern wind power facilities are developing towards the trend of gigantism and dense layout. This is not only reflected in the expansion of the scale of single wind turbines, but also the planning of the entire wind farm tends to be compact, striving to deploy more equipment within a limited geographical area. The resulting problem is the impact of wake effects on the overall performance of the wind farm. And accurately describing the meandering of the yaw fan wake is crucial. In past studies, most yaw analysis models used superposition models to consider the interaction of wind turbine wakes, which would bring certain errors. The present invention proposes a new yaw fan wake calculation model method, based on the CFD actuator line model and an improved Jensen wake model, for analyzing the yaw of tandem wind turbines. When the wind turbine yaws or moves, the center of its wake will inevitably shift. As the wake develops, it will cause large-scale expansion and meandering of the wake, which will cause great difficulties in wake prediction. Accurate tracking of the wake center is crucial. Therefore, the present invention proposes a search zero-point method for the problems existing in the application of the existing wake center tracking method to the wake model. It solves the problems of the existing wake center tracking method and makes it applicable to the hybrid wake model. This not only realizes wake center tracking but also ensures the velocity continuity between hybrid models. The aim is to accurately calculate the impact of the wind turbine wake considering the yaw angle on the wind farm. Summary of the Invention

[0003] The purpose of the present invention is to provide a new yaw fan wake calculation model method, introducing a search zero-point method, which solves the problems of insufficient accuracy and applicability of the existing wake model in wake center tracking. In addition, the present invention also uses the tracer particle method to calculate the meandering of the yaw wake, effectively solving the problem of predicting the wind turbine wake considering the yaw factor. The specific implementation is as follows: Use the CFD method based on the actuator line theory to simulate and calculate the near wake field of the wind turbine, then use the search zero-point method to determine the wake center, and then use the improved Jensen model to simulate and calculate the far wake field of the wind turbine. In the coupling area, the exponential factor is used to neutralize and smooth the calculation of CFD and Jensen. Compared with the medium-precision wake model, this coupling method has higher calculation accuracy; compared with the high-precision model, this coupling method has lower calculation cost and is suitable for preliminary estimation calculations in engineering.

[0004] To achieve the above object, the present invention provides a new method for calculating the wake of a yaw wind turbine, which includes the following steps:

[0005] S1. When performing numerical simulation of the wind turbine wake effect, divide the computational domain into three regions, including the near-wake region C1 of the wind turbine, the coupling region C2, and the far-wake region C3 of the wind turbine. Set the time step to T, and set the initial conditions for the velocity fields of the three regions: V C1 , V C2 and V C3 . Initialize the velocity vector as (V wind , 0, 0), where V wind represents the wind speed of the initial uniform wind field.

[0006] S2. During the numerical simulation of the near-wake region C1 of the wind turbine, establish a virtual actuator line model of the wind turbine based on the actuator line theory. Based on the Reynolds-averaged NS equations, introduce a source term to replace the actual wind turbine blades to obtain the velocity field V C11 in the near-wake region of the wind turbine.

[0007] S3. When the wind turbine yaws or moves, it affects the wake. Use the search zero-point method to determine the wake center.

[0008] S4. Use the tracer particle method to calculate the yaw wake meandering.

[0009] S5. During the numerical simulation of the coupling region C2 of the wind turbine, add an exponential relaxation factor to perform coupling calculations on the velocity fields V C11 and V C3 to obtain the velocity field V C2m in the coupling region, where m = 1, 2, 3,....

[0010] S6. During the numerical simulation of the far-wake region C3 of the wind turbine, based on the existing improved Jensen model, combine the velocity field V C2m calculated in step S5 with the wake center P n obtained in S4 as the inlet of the improved Jensen model to obtain the boundary conditions, and then calculate the far-wake velocity field V C3m of the wind turbine through the improved Jensen model.

[0011] S7. Use the new V C31 to repeat step S5 and output the new velocity field V C2(m+1) in the new coupling region.

[0012] S8. Repeat S6 and S7 until the final required far-wake velocity field V C3(m+1) of the wind turbine is calculated.

[0013] Preferably, in step S2, the velocity field V of the near-wake region of the wind turbine is calculated as follows: C11 The process is as follows:

[0014] Set the length of the computational domain of the near-wake region C1 of the wind turbine to x C1 , divide the computational domain of the near-wake region C1 of the wind turbine into grids, and solve the initial velocity field V of the near-wake region of the wind turbine C1 using the Reynolds-averaged NS equations, and add a source term in the Reynolds-averaged NS equations using the actuator line model of the wind turbine to replace the real wind turbine blades. The formula is as follows:

[0015]

[0016]

[0017] In the formula, u is the velocity vector of the near-wake computational grid, u j is the velocity component in the j direction, x j is the axial coordinate axis component of the wind field, ρ is the fluid density, u i is the velocity component in the i direction, t is the wind field calculation time, p is the pressure, τ ij and are the viscous stress and the turbulent stress, g i is the gravitational acceleration, f σi is the surface tension, f UALM is the source term for modeling the blade effect on the fluid, that is, the vector sum of the lift and drag. i = 1, 2, 3 represent the x, y, and z coordinate axes in the Cartesian coordinate system respectively, and j = 1, 2, 3 represent the x, y, and z coordinate axes in the Cartesian coordinate system respectively;

[0018] In the near-wake field of the wind turbine, using the actuator line model of the wind turbine, each blade of the wind turbine is discretized into a series of spanwise sections with fixed profiles, chord lengths, and angles of attack, and each spanwise section is replaced by an actuator point that generates a body force. The body force on these actuator points is represented by f k . The aerodynamic force generated by each blade segment at (x, y, z, t) is given by the following formula:

[0019]

[0020] In the above formula, F(x, y, z, t) is the vector sum of the lift and drag of each aerodynamic point, x, y, z, and t are the three-dimensional coordinate points and the current wind field calculation time respectively, n is the actuator point index, N is the total number of actuator point segments, d n is the distance between (x, y, z) and the position of the nth actuator point, and ε is a constant that determines the width of the projection area;

[0021] The formula for the velocity field of the near-wake region of the wind turbine is obtained as follows:

[0022] V C11 = (v1, v2, v3);

[0023] When the time step is reached, the near wake flow field calculation is completed, where v1, v2, and v3 represent the velocities in the x, y, and z coordinate axes directions in the Cartesian coordinate system respectively.

[0024] Preferably, in step S3, the search for the wake center using the zero search method is as follows:

[0025] S31. Use the CFD method to obtain the wake center P0 at t = x C1 / V wind ;

[0026] S32. When t > x C1 / V wind , use P0 as the geometric center to find the zero point of the wake tangential velocity within a square region with a length and width of 0.3R. The zero point of Vz determines the Y-axis coordinate of the wake center, and the zero point of Vy determines the Z-axis coordinate of the wake center, obtaining a new wake center P1;

[0027] S33. In the next time step, use P1 as the new wake center and repeat step S32 to obtain the next new wake center P2;

[0028] S34. By iteratively repeating steps S31 to S33, track the zero points of the tangential velocity within the wake at each time step to obtain P n , and then regard P n as the wake center of the Jensen inlet boundary condition; obtain the wake center cylindrical coordinates (x g , y g , z g ).

[0029] Preferably, in step S4, the calculation of the yawed wake meandering specifically includes:

[0030] Design a quasi-Lagrangian method to describe the large eddy turbulence mechanism. Among them, the simulated eddies are larger than two rotor diameters. Set v and w as the lateral and vertical velocities respectively, and define the characteristic velocities as follows:

[0031]

[0032] Among them, (x b , y b , z b ) is the inertial coordinate system, A f is the wake cross-sectional area, where the central coordinates of the wake cross-section are (y b , z b ), vc is the lateral characteristic velocity, w c is the vertical characteristic velocity;

[0033] The lateral and vertical wake displacements are described as:

[0034]

[0035] where y g is the lateral wake displacement, U is the incoming flow wind speed, z g is the vertical wake displacement, according to the quasi-Lagrangian method T m = L / U, T m is the time interval, L is the wake length under consideration, t i is the time when the wake is released. In addition, the initial condition of the wake release time t i is given as follows:

[0036]

[0037] Preferably, in step S4, based on the theory of the quasi-Lagrangian method, the specific solution of the wake meandering process is as follows:

[0038] The time t q = t0 + qΔt; q = 0, 1, 2, 3,..., R; where t0 is the initial time, Δt is the time step, and for each time step, the following solution process is as follows:

[0039] Determine the initial conditions, determine the lateral and vertical wake displacements, and determine the wake center Repeat determining the lateral and vertical wake displacements and the wake center until k + 1 = R, where the wake center is:

[0040]

[0041] Preferably, in step S5, the calculation process of the coupled region velocity field V C2m is as follows:

[0042] The velocity at each point in the coupled region is calculated by coupling two velocity fields V C2 and V C3 to make the velocity distribution in the coupled region smooth. The calculation formula is as follows:

[0043] V C2m = V C11 ω + V C3 (1 - ω);

[0044] In the formula, ω is the weight factor:

[0045]

[0046] In the formula, δ is the dimensionless length or relative length of the coupling region C2, and the formula is as follows:

[0047]

[0048] In the above formula, H is the length of the computational domain from the starting point to the ending point of the wake at the current time, x C1 is the length of the computational domain of the near-wake region C1, and x C2 is the length of the computational domain of the coupling region C2.

[0049] Preferably, in step S6, the velocity V C3m in the far-wake region is calculated, where V C3m = Ar 2 + Br + C, where A, B, and C are constants respectively, and satisfy the following conditions:

[0050] When the radial radius expands from r to r x , the wake velocity recovers to the natural wind speed, and the wake velocity at the highest point and the lowest point are:

[0051]

[0052] In the formula, r x is the cross-sectional radius of the far-wake region, β is the wind shear coefficient, V0 is the velocity of the wake center at the entrance of the far-wake region, and h is the height of the wake center at the entrance of the far-wake region;

[0053] If the wake radius of the improved Jensen wake model is the same as that of the Jensen wake model and the mass fluxes are equal, then:

[0054]

[0055] In the formula, u is the velocity in the wake region of the original Jensen model;

[0056] Preferably, the velocity solution in the far-wake region is carried out in two steps:

[0057] Calculation and prediction step of the original Jensen model:

[0058] k = 0.5 / ln(h / z0)

[0059]

[0060] u = u0[1 - 2a / (1 + kx c3 / r0) 2 ;

[0061] In the formula, k represents the wake expansion coefficient, z0 represents the ground roughness, C Trepresents the thrust coefficient, a represents the axial induction factor, x C3 represents the length of the far wake region calculation domain, u0 represents the incoming flow wind speed, and r0 represents the initial wake radius;

[0062] Correction step:

[0063]

[0064] Therefore, the present invention adopts the above-mentioned new yaw wind turbine wake calculation model method, and the technical effects are as follows:

[0065] (1) In the present invention, a zero-point search method is proposed, which solves the problems of low accuracy and inapplicability of the existing wake center tracking method in the hybrid wake model. The zero-point search method ensures the continuity of the velocity between the two wake models and successfully simulates the wake deflection of the yaw wind turbine.

[0066] (2) In the present invention, the meandering of the yaw wake is considered. Based on the dynamic wake meandering model, the wake meandering is simulated and calculated by using the tracer particle method, which solves the problem that the existing wake model does not consider the meandering of the yaw wake.

[0067] (3) In the present invention, the problem of velocity backflow can be avoided even when CFD uses a small calculation domain through the method of numerical relaxation, ensuring the stability of the velocity-pressure coupling calculation in the calculation domain. Description of the drawings

[0068] Figure 1 It is a flow chart of a new yaw wind turbine wake calculation model method of the present invention. Specific implementation manners

[0070] The technical solutions of the present invention will be further described below with reference to the drawings and embodiments.

[0071] Unless otherwise defined, the technical terms or scientific terms used in the present invention should have the ordinary meaning understood by those of ordinary skill in the field to which the present invention belongs.

[0072] Embodiment 1

[0073] A new yaw wind turbine wake calculation model method includes the following steps:

[0074] S1. When performing the numerical simulation of the wind turbine wake effect, the calculation domain is divided into three regions, including the wind turbine near wake region C1, the coupling region C2, and the wind turbine far wake region C3. The time step is set to T, and the initial conditions of the velocity fields of the three regions are set: V C1 、V C2 and V C3 , and the velocity vector is initialized to (Vwind , 0, 0), where V wind represents the wind speed of the initial uniform wind field;

[0075] S2. During the numerical simulation of the near-wake region C1 of the wind turbine, an actuator line model of the wind turbine is established based on the actuator line theory. Based on the Reynolds-averaged NS equations, a source term is introduced to replace the actual wind turbine blades to obtain the velocity field V of the near-wake region of the wind turbine C11 ;

[0076] Set the computational domain length of the near-wake region C1 of the wind turbine to x C1 , divide the computational domain of the near-wake region C1 of the wind turbine into grids, and use the Reynolds-averaged NS equations to solve the velocity field V of the initial near-wake region of the wind turbine C1 Add a source term in the Reynolds-averaged NS equations using the actuator line model of the wind turbine to replace the real wind turbine blades. The formula is as follows:

[0077]

[0078] In the formula, u is the velocity vector of the near-wake computational grid, u j is the velocity component in the j direction, x j is the axial coordinate axis component of the wind field, ρ is the fluid density, u i is the velocity component in the i direction, t is the wind field calculation time, p is the pressure, τ ij and are the viscous stress and turbulent stress, g i is the gravitational acceleration, f σi is the surface tension, f UALM is the source term for modeling the blade effect on the fluid, that is, the vector sum of the lift and drag. i = 1, 2, 3 represent the x, y, z coordinate axes in the Cartesian coordinate system respectively, and j = 1, 2, 3 represent the x, y, z coordinate axes in the Cartesian coordinate system respectively;

[0079] In the near-wake flow field of the wind turbine, using the actuator line model of the wind turbine, each blade of the wind turbine is discretized into a series of spanwise sections with fixed profiles, chord lengths, and angles of attack. Each spanwise section is replaced by an actuator point that generates a body force. The body force on these actuator points is represented by f k . The aerodynamic force generated by each blade segment at (x, y, z, t) is given by the following formula:

[0080]

[0081] In the above formula, F(x, y, z, t) is the vector sum of the lift and drag of each aerodynamic point. x, y, z, t are the three-dimensional coordinate points and the current wind field calculation time respectively. n is the actuator point index, N is the total number of actuator point segments, d nis the distance between (x, y, z) and the nth actuator point position, and ε is a constant determining the width of the projection area;

[0082] The velocity field formula in the near wake region of the wind turbine is obtained as follows:

[0083] V C11 = (v1, v2, v3);

[0084] When the time step is reached, the calculation of the near wake field domain is completed, and v1, v2, and v3 respectively represent the velocities in the x, y, and z coordinate axis directions in the Cartesian coordinate system.

[0085] S3. The influence on the wake caused by the yaw or movement of the wind turbine. The search zero point method is used to determine the wake center, specifically as follows:

[0086] S31. Use the CFD method to obtain the wake center P0 at t = x C1 / V wind ;

[0087] S32. When t > x C1 / V wind , use P0 as the geometric center to find the zero point of the tangential velocity of the wake within a square area with a length and width of 0.3R. The zero point of Vz determines the Y-axis coordinate of the wake center, and the zero point of Vy determines the Z-axis coordinate of the wake center, obtaining a new wake center P1;

[0088] S33. In the next time step, use P1 as the new wake center and repeat step S32 to obtain the next new wake center P2;

[0089] S34. By iteratively repeating steps S31 to S33, track the zero point of the tangential velocity within the wake at each time step, thereby obtaining P n , and then regard P n as the wake center of the Jensen inlet boundary condition; obtain the wake center cylindrical coordinates (x g , y g , z g ).

[0090] S4. Use the tracer particle method to calculate the yaw wake meandering, specifically including:

[0091] Design a quasi-Lagrangian method to describe the large eddy turbulence mechanism. Among them, the simulated eddies are larger than two rotor diameters. Set v and w as the transverse and vertical velocities respectively, and define the characteristic velocities as follows:

[0092]

[0093] Among them, (x b , y b, z b ) is an inertial coordinate system, A f is the wake cross-sectional area, where the origin is (y b , z b ) has a diameter D w , v c is the transverse characteristic velocity, w c is the vertical characteristic velocity;

[0094] The transverse and vertical wake displacements are described as:

[0095]

[0096] where, y g is the transverse wake displacement, U is the incoming flow wind speed, z g is the vertical wake displacement, according to the quasi-Lagrangian method T m = L / U, T m is the time interval, L is the wake length under consideration, t i is the time when the wake is released, and in addition, the initial condition of the wake release time t i is given as follows:

[0097]

[0098]

[0099] Based on the theory of the quasi-Lagrangian method, the specific solution of the wake meandering process is as follows:

[0100] Represent the velocity field at the Jensen inlet boundary in cylindrical coordinates as V C = (V x , V y , V z );

[0101] The time t q = t0 + qΔt; q = 0, 1, 2, 3,..., R; where t0 is the initial time, Δt is the time step, and for each time step, there is the following solution process:

[0102] Determine the initial conditions, determine the transverse and vertical wake displacements, and determine the wake center Repeat determining the transverse and vertical wake displacements and the wake center until k + 1 = R, where the wake center is:

[0103]

[0104] S5. During the numerical simulation of the wind turbine coupling region C2, an exponential relaxation factor is added to the velocity fields V C11 and V C3Perform coupling calculations to obtain the velocity field V in the coupling region C2m , m = 1, 2, 3,....;

[0105] The velocity at each point in the coupling region is calculated by coupling two velocity fields V C2 and V C3 to smooth the velocity distribution in the coupling region. The calculation formula is as follows:

[0106] V C2m = V C11 ω + V C3 (1 - ω);

[0107] In the formula, ω is the weight factor:

[0108]

[0109] In the formula, δ is the dimensionless length or relative length of the coupling region C2, and the formula is as follows:

[0110]

[0111] In the above formula, H is the length of the computational domain from the starting point to the ending point of the wake at the current time, x C1 is the length of the computational domain of the near-wake region C1, and x C2 is the length of the computational domain of the coupling region C2.

[0112] S6. During the numerical simulation of the near-wake region C3 of the wind turbine, based on the existing improved Jensen model, the velocity field V C2m calculated in step S5 is combined with the wake center P n obtained in S4 as the inlet of the improved Jensen model to obtain the boundary conditions, and then the far-wake velocity field V C3m of the wind turbine is calculated through the improved Jensen model;

[0113] Calculate the velocity V C3m in the far-wake region, where V C3m = Ar 2 + Br + C, and A, B, and C are constants respectively and satisfy the following conditions:

[0114] (1) When the radial radius expands from r to r x , the wake velocity recovers to the natural wind speed, and the highest and lowest wake velocities are:

[0115]

[0116] In the formula, r xwhere \(r_0\) is the cross-sectional radius of the far wake region, \(\beta\) is the wind shear coefficient, \(V_0\) is the velocity at the center of the wake at the entrance of the far wake region, and \(h\) is the height of the center of the wake at the entrance of the far wake region;

[0117] (2) If the wake radius of the improved Jensen wake model is the same as that of the Jensen wake model and the mass fluxes are equal, then:

[0118]

[0119] where \(u\) is the velocity in the wake region of the original Jensen model;

[0120] The velocity solution in the far wake region is carried out in two steps:

[0121] (1) Prediction step by the original Jensen model:

[0122] \(k = 0.5 / \ln(h / z_0)\)

[0123]

[0124] where \(k\) represents the wake expansion coefficient, \(z_0\) represents the ground roughness, \(C\) T represents the thrust coefficient, \(a\) represents the axial induction factor, \(x\) C3 represents the length of the calculation domain in the far wake region, \(u_0\) represents the incoming flow velocity, and \(r_0\) represents the initial wake radius;

[0125] (2) Correction step:

[0126]

[0127] S7. Use the new \(V\) C31 Repeat step S5 and output the new velocity field \(V\) in the coupling region C2(m+1) ;

[0128] S8. Repeat S6 and S7 until the velocity field \(V\) of the far wake region of the wind turbine required finally is calculated C3(m+1) .

[0129] Therefore, the present invention adopts the above-mentioned new yaw wind turbine wake calculation model method, aiming to improve the accuracy and efficiency of calculation. By combining the CFD actuator line model and the improved Jensen wake model, and introducing a search zero-point method, the problems of insufficient accuracy and applicability of the existing wake model in wake center tracking are solved. In addition, the present invention also uses the tracer particle method to calculate the meandering of the yaw wake, effectively solving the problem of predicting the wind turbine wake considering the yaw factor. The present invention uses the CFD method based on the actuator line theory to simulate and calculate the near wake field of the wind turbine, then uses the search zero-point method to determine the wake center, and then uses the improved Jensen model to simulate and calculate the far wake field of the wind turbine. In the coupling area, the exponential factor is used to neutralize the CFD and Jensen for smooth calculation. Compared with the medium-precision wake model, this coupling method has higher calculation accuracy; compared with the high-precision model, this coupling method has lower calculation cost.

[0130] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit them. Although the present invention has been described in detail with reference to the preferred embodiments, those of ordinary skill in the art should understand that they can still modify or equivalently replace the technical solutions of the present invention, and these modifications or equivalent replacements cannot make the modified technical solutions deviate from the spirit and scope of the technical solutions of the present invention.

Claims

1. A novel yaw wind turbine wake calculation model method, characterized in that: The following steps are involved: S1. When performing numerical simulation of the wind turbine wake effect, the calculation domain is divided into three regions, including the wind turbine near wake region C1, the coupling region C2, and the wind turbine far wake region C3. The time step is set to T, and the initial conditions are set to the velocity fields of the three regions: V C1 、V C2 and V C3 , the velocity vector is initialized as (V wind ,0,0), where V wind represents the wind speed of the initial uniform wind field; S2. In the process of numerical simulation of the near-wake area C1 of the wind turbine, the actuating line model of the wind turbine is established based on the actuating line theory. Based on the Reynolds averaged NS equation, the source term is introduced to replace the actual wind turbine blades to obtain the velocity field V in the near-wake area of ​​the wind turbine. C11 ; S3. When the wind turbine yaws or moves, it affects the wake, and the search zero method is used to determine the wake center; S4, Calculate the yaw wake meander using the tracer particle method; S5. In the numerical simulation of the wind turbine coupling region C2, an exponential relaxation factor was added to the velocity field V C11 and V C3 Perform coupling calculation to obtain the velocity field V in the coupling zone C2m , m=1,2,3,....; S6. In the process of numerical simulation of the near-wake area C3 of the wind turbine, based on the existing improved Jensen model, the velocity field V calculated in step S5 is C2m Combined with the wake center P obtained in S4 n As the entry of the improved Jensen model, the boundary conditions are obtained, and then the wind turbine far wake velocity field V is calculated by the improved Jensen model. C3m ; S7, use the new V C31 Repeat step S5 to output the new coupling zone velocity field V C2(m+1) ; S8, repeat S6 and S7 until the final required velocity field V in the far wake area of ​​the wind turbine is calculated C3(m+1) .

2. A novel yaw wind turbine wake calculation model method according to claim 1, characterized in that: In step S2, the velocity field V in the wind turbine wake region is calculated. C11 The process is as follows: Set the computational domain length of the wind turbine wake area C1 to x C1 , divide the computational domain of the wind turbine wake area C1 into grids, and calculate the initial wind turbine wake area velocity field V C1 The Reynolds averaged NS equation is used to solve the problem. The fan actuation line model is used to add source terms to the Reynolds averaged NS equation to replace the real fan blades. The formula is as follows: Where u is the velocity vector of the near wake calculation grid, u j is the velocity component in the j direction, x j is the axial coordinate axis component of the wind field, ρ is the fluid density, u i is the velocity component in the i direction, t is the wind field calculation time, p is the pressure, τ ij and are the viscous stress and turbulent stress, g i is the acceleration due to gravity, f σi is the surface tension, f UALM is the source term for modeling the blade effect on the fluid, i.e., the vector sum of lift and drag, i=1,2,3 respectively represent the directions of the x, y, and z axes in the Cartesian coordinate system, j=1,2,3 respectively represent the directions of the x, y, and z axes in the Cartesian coordinate system; In the near-wake flow field of a wind turbine, the wind turbine actuation line model is used to discretize each blade of the wind turbine into a series of spanwise sections with fixed profiles, chord lengths and angles of attack. Each spanwise section is replaced by an actuation point that generates body force. The body force on these actuation points is expressed by f k The aerodynamic formula generated by each blade segment at (x, y, z, t) is as follows: In the above formula, F(x,y,z,t) is the vector sum of lift and drag of each aerodynamic point, x, y, z, t are the three-dimensional coordinate point and the current wind field calculation time, n is the actuation point index, N is the total number of actuation point segments, d n is the distance between (x, y, z) and the position of the nth actuator point, and ε is a constant that determines the width of the projection area; The formula for the velocity field near the wind turbine wake area is as follows: V C11 =(v1,v2,v3); When the time step When , the calculation of the near wake field is completed, and v1, v2, and v3 represent the velocities in the x, y, and z axes in the Cartesian coordinate system, respectively.

3. A novel yaw wind turbine wake calculation model method according to claim 1, characterized in that: In step S3, the wake center is determined using the zero-point search method, as follows: S31. Using CFD method to obtain t=x C1 / V wind The wake center at P0; S32, at t>x C1 / V wind When P0 is used as the geometric center, the zero point of the wake tangential velocity is found in a square area with a length and width of 0.3R. The zero point of Vz determines the Y-axis coordinate of the wake center, and the zero point of Vy determines the Z-axis coordinate of the wake center, and the new wake center P1 is obtained; S33, in the next time step, use P1 as the new wake center, repeat step S32, and obtain the next new wake center P2; S34, by iteratively repeating steps S31 to S33, tracking the tangential velocity zero point in the wake at each time step, thereby obtaining P n , then P n The wake center is regarded as the Jensen inlet boundary condition; the wake center cylindrical coordinates (x g ,y g ,z g ).

4. A novel yaw wind turbine wake calculation model method according to claim 1, characterized in that: In step S4, the yaw wake meander is calculated, which specifically includes: A quasi-Lagrangian method is designed to describe the large eddy turbulence mechanism, where the simulated eddy is larger than two rotor diameters, v and w are set as the lateral and vertical velocities respectively, and the characteristic velocity is defined as follows: Among them, (x b ,y b ,z b ) is the inertial coordinate system, A f is the wake cross-sectional area, where the center coordinate of the wake cross-sectional area is (y b ,z b ), v c is the lateral characteristic velocity, w c is the vertical characteristic velocity; The lateral and vertical wake displacements are described by: Among them, y g is the lateral wake displacement, U is the incoming wind speed, z g is the vertical wake displacement, according to the quasi-Lagrangian method T m =L / U,T m is the time interval, L is the wake length considered, t i is the time of wake release, and the time of wake release t is also given. i The initial conditions are as follows: y g (t i )=0; z g (t i )=0; 5. A novel yaw wind turbine wake calculation model method according to claim 4, characterized in that: In step S4, based on the theory of the quasi-Lagrangian method, the wake meandering process is specifically solved as follows: Time t q =t0+qΔt; q=0,1,2,3,…,R; where t0 is the initial time, Δt is the time step, and the solution process for each time step is as follows: Determine initial conditions, determine lateral and vertical wake displacements, and determine the wake center Repeat the determination of the lateral and vertical wake displacements and the wake center until k+1=R, where the wake center is:

6. A novel yaw wind turbine wake calculation model method according to claim 1, characterized in that: In step S5, the coupling zone velocity field V is calculated C2m The calculation process is as follows: The velocity of each point in the coupling zone is calculated using two velocity fields V C2 and V C3 Carry out coupling calculation to make the velocity distribution in the coupling zone smooth. The calculation formula is as follows: V C2m =V C11 ω+V C3 (1-ω); Where ω is the weight factor: In the formula, δ is the dimensionless length or relative length of the coupling region C2, and the formula is as follows: In the above formula, H is the length of the calculation domain from the start point to the end point of the tail flow at the current time, x C1 is the length of the computational domain near the wake region C1, x C2 is the calculation domain length of the coupling region C2.

7. A novel yaw wind turbine wake calculation model method according to claim 1, characterized in that: In step S6, the velocity V in the far wake region is calculated. C3m , where V C3m =Ar 2 +Br+C, A, B, C are constants respectively, and the following conditions are met: When the radial radius expands from r to r x When , the wake velocity returns to the natural wind speed, then the highest wake velocity and the lowest wake velocity are: In the formula, r x is the cross-sectional radius of the far wake region, β is the wind shear coefficient, V0 is the velocity of the wake center at the entrance of the far wake region, and h is the height of the wake center at the entrance of the far wake region; The improved Jensen wake model and the Jensen wake model have the same wake radius and the same mass flux, so: Where u is the wake velocity of the original Jensen model.

8. A novel yaw wind turbine wake calculation model method according to claim 7, characterized in that: The velocity solution in the far wake region is carried out in two steps: The original Jensen model calculates the prediction step: k=0.5 / ln(h / z0) u=u0[1-2a / (1+kx c3 / r0) 2 ]; In the formula, k represents the wake expansion coefficient, z0 represents the ground roughness, C T represents the thrust coefficient, a represents the axial induction factor, x C3 represents the length of the calculation domain of the far wake area, u0 represents the incoming wind speed, and r0 represents the initial wake radius; Correction step:

Images