Rotary fan blade water drop collision coefficient determination method based on Lagrangian method

Abstract

The invention provides a rotating fan blade water drop collision coefficient determination method based on a Lagrangian method, and the method comprises the following steps: S1, obtaining environmental parameters and a fan rotating speed in an actual working condition, the environmental parameters comprising an incoming flow wind speed, the size of liquid drops in air, and atmospheric pressure; s2, simulating fan rotation by adopting a motion reference system method based on the environmental parameters and the fan rotation speed in the actual working condition, and establishing an air flow field model; s3, on the basis of the air flow field model, turbulent kinetic energy and specific dissipation rate are determined, and derivatives of the turbulent kinetic energy and the specific dissipation rate to time are determined; s4, constructing an air water drop phase control equation; s5, the turbulent kinetic energy, the derivative of the specific dissipation rate to time and the water drop phase control equation are input into a DPM model in FLUENT software to determine the distance between the adjacent water drops at the release position and the distance between the adjacent water drops hitting the fan blade; and S6, determining a water drop collision coefficient beta.

Description

Technical Field

[0001] The present invention relates to the field of wind turbines, and particularly to a method for determining the water droplet collision coefficient of a rotating wind turbine blade based on the Lagrangian method. Background Art

[0002] Under icing meteorological conditions, there are usually a large number of supercooled water droplets in the air. Since the water droplets lack condensation nuclei in the air, it is difficult to freeze even when the temperature is below 0°C. These supercooled water droplets are usually extremely unstable and are prone to phase change under disturbance. When the wind turbine is operating normally, the supercooled water droplets suspended in the air hit the blade surface and heterogeneous nucleation occurs on the cold surface, thereby triggering blade icing, and the icing of the wind turbine (abbreviation: wind turbine or wind power generator) blade will reduce the operating stability of the wind turbine.

[0003] The water droplet collision coefficient reflects the occurrence and development process of blade icing. In the prior art, the water droplet collision coefficient is often determined by the velocities before and after the collision. However, in actual working conditions, it is difficult to accurately determine the velocities of the water droplets before and after the collision. Moreover, since the rotation process of the wind turbine blade will also affect the water droplet collision, the water droplet collision coefficient of the wind turbine blade cannot be accurately determined.

[0004] Therefore, in order to solve the above technical problems, it is urgent to propose a new technical means. Summary of the Invention

[0005] In view of this, the purpose of the present invention is to provide a method for determining the water droplet collision coefficient of a rotating wind turbine blade based on the Lagrangian method. By simulating and constructing the air flow field of the wind turbine blade through the environmental parameters of the actual working condition environment and the rotational speed of the wind turbine, fully considering the air flow disturbance effect caused by the blade rotation, and calculating the water droplet collision rate based on the Lagrangian method, it can effectively ensure the accuracy of the final calculation result and provide accurate data support for subsequent icing warning and de-icing measure formulation.

[0006] A method for determining the water droplet collision coefficient of a rotating wind turbine blade based on the Lagrangian method provided by the present invention includes the following steps:

[0007] S1. Obtain the environmental parameters and the rotational speed of the wind turbine in the actual working condition, where the environmental parameters include the incoming flow wind speed, the size of the liquid droplets in the air, and the atmospheric pressure;

[0008] S2. Simulate the rotation of the wind turbine by using the moving reference frame method based on the environmental parameters and the rotational speed of the wind turbine in the actual working condition, and establish an air flow field model;

[0009] S3. Determine the turbulent kinetic energy and the specific dissipation rate based on the air flow field model, and determine the derivatives of the turbulent kinetic energy and the specific dissipation rate with respect to time;

[0010] S4. Analyze the forces acting on the droplets in the gas-liquid two-phase flow field using the Lagrangian method, and construct the governing equation for the water droplet phase in the air;

[0011] S5. Input the time derivatives of the turbulent kinetic energy and the specific dissipation rate, as well as the governing equation for the water droplet phase, into the DPM model in the FLUENT software to determine the distance between adjacent droplets at the release position and the distance between adjacent droplets hitting the fan blade;

[0012] S6. Determine the droplet collision coefficient β:

[0013]

[0014] where: dY represents the distance between adjacent droplets at the release position, dL represents the distance between adjacent droplets hitting the fan blade, β1 is the droplet collision coefficient on the two-dimensional plane of the fan blade airfoil, and z1 and z2 are the distances of the adjacent droplet collision blade positions from the origin of the z-axis coordinate in the established three-dimensional coordinate system.

[0015] Furthermore, the air flow field model is:

[0016]

[0017] where: represents the air velocity; u, v, and w respectively represent the components of the air velocity in the x, y, and z directions, p represents the atmospheric pressure; υ represents the air dynamic viscosity.

[0018] Furthermore, determining the turbulent kinetic energy and the specific dissipation rate based on the air flow field model, and determining the time derivatives of the turbulent kinetic energy and the specific dissipation rate specifically include:

[0019] Construct the mass equation of the three-dimensional air flow field in the Cartesian coordinate system:

[0020] where i, j = x, y, z; ρ is the air density, μ is the kinematic viscosity coefficient, represents the turbulent stress;

[0021] where: P t is the pressure caused by the pulsating velocity; where: k is the turbulent kinetic energy per unit mass of the fluid;

[0022] δ ij is the Dirac function, and the expression is:

[0023]

[0024] μt represents the turbulent viscosity coefficient in the k-ω equation at low Reynolds numbers:

[0025]

[0026] In the k-ε equation at high Reynolds numbers:

[0027]

[0028] Determine the derivatives of the turbulent kinetic energy k and the specific dissipation rate ω with respect to time t:

[0029] where P k represents the pressure caused by turbulence; ω is the specific dissipation rate; k is the turbulent kinetic energy; β * , σ k are closure constants; μ t is the turbulent viscosity coefficient.

[0030]

[0031] where: F1 and F2 are auxiliary functions respectively, where:

[0032]

[0033] Furthermore, constructing the control equation for the water droplet phase in the air specifically includes:

[0034]

[0035] In the rectangular coordinate system, rewrite it in matrix form:

[0036]

[0037] where: K is the air-water droplet exchange coefficient, F x is the additional force caused by rotation;

[0038]

[0039] C D is the drag coefficient, R e represents the relative Reynolds number;

[0040]

[0041] The beneficial effects of the present invention: Through the present invention, the air flow field of the fan blade is simulated and constructed based on the environmental parameters of the actual working condition environment and the rotational speed of the fan, fully considering the air flow disturbance effect caused by the rotation of the blade, and calculating the water droplet collision rate based on the Lagrangian method, which can effectively ensure the accuracy of the final calculation result and provide accurate data support for subsequent icing warning and de-icing measure formulation. Description of the Drawings

[0042] The present invention will be further described below in conjunction with the drawings and embodiments:

[0043] Figure 1 It is a schematic flow diagram of the present invention.

[0044] Figure 2 It is the blade grid division.

[0045] Figure 3 It is a schematic diagram of the stationary and rotating coordinate systems.

[0046] Figure 4 It is a schematic diagram of the blade spanwise direction.

[0047] Figure 5 It is the airfoil diagram of the blade cross-section. Detailed Description of the Invention

[0048] The following further details the present invention:

[0049] A method for determining the water droplet collision coefficient of a rotating fan blade based on the Lagrangian method provided by the present invention includes the following steps:

[0050] S1. Obtain the environmental parameters and the fan speed in the actual working condition. Among them, the environmental parameters include the incoming flow wind speed, the liquid droplet size in the air, and the atmospheric pressure;

[0051] S2. Simulate the rotation of the fan using the moving reference frame method based on the environmental parameters and the fan speed in the actual working condition, and establish an air flow field model; among them, the moving reference frame method (abbreviated as the MRF method) is used to simulate the rotation of the fan. First, determine the shape of the fan blade, construct a simulation model, and perform grid division based on the finite element idea. The grid is a microelement control volume, and then the MRF method is used to simulate the rotation;

[0052] S3. Determine the turbulent kinetic energy and the specific dissipation rate based on the air flow field model, and determine the derivatives of the turbulent kinetic energy and the specific dissipation rate with respect to time;

[0053] S4. Construct the control equation for the water droplet phase in the air;

[0054] S5. Input the derivatives of the turbulent kinetic energy and the specific dissipation rate with respect to time and the control equation for the water droplet phase into the DPM model in the FLUENT software to determine the distance between adjacent water droplets at the release position and the distance between adjacent water droplets hitting the fan blade;

[0055] S6. In the three-dimensional fan blade model shown in Figure 2 , according to the required blade spanwise Z coordinate, determine the two-dimensional airfoil of the blade cross-section and the water droplet collision situation at the cross-section, extract the spatial calculation range and boundary conditions, refer toFigure 4 Based on the blade model, corresponding blade cross-sections can be extracted on different z-planes, thereby obtaining the geometric shape of the two-dimensional airfoil, as Figure 5 shown.

[0056] Determine the water droplet collision coefficient β:

[0057]

[0058] where: dY represents the distance between adjacent water droplets at the release position, dL represents the distance between adjacent water droplets impinging on the wind turbine blade, β1 is the water droplet collision coefficient on the two-dimensional plane of the wind turbine blade airfoil, and z1 and z2 are the distances of the positions where adjacent water droplets collide with the blade relative to the origin of the z-axis coordinate in the established three-dimensional coordinate system, as shown in the appendix Figure 4 By simulating and constructing the air flow field of the wind turbine blade through the environmental parameters of the actual working condition environment and the rotational speed of the wind turbine, fully considering the air flow disturbance effect caused by the blade rotation, and calculating the water droplet collision rate based on the Lagrangian method, the accuracy of the final calculation result can be effectively ensured, providing accurate data support for subsequent icing warning and de-icing measure formulation.

[0059] Simulate the rotation of the wind turbine blade based on the MRF method. This method uses a reference coordinate system with acceleration to characterize the rotational flow. When processing the computational domain of the numerical simulation, it is divided into two parts: a rotating domain and a stationary domain. Among them, the area far from the outer flow field of the wind turbine blade remains the stationary area, and the control equations of its external flow field are solved in the inertial reference coordinate system; while a spanwise acceleration is introduced in the near-wall area of the blade, and the flow field is solved according to the control equations in the rotating reference system. Data transfer between the two areas is carried out through the coupling of the interface. Define a rotating reference system that follows the blade rotation, with the hub of the wind turbine as the coordinate origin, the rotation axis coinciding with the wind turbine rotation axis, and the rotational speed being the same as the wind turbine rotational speed. In this way, the flow field is transformed into a steady state, and then the control equations in the relative coordinate system are solved. In addition, on the domain interface, the flow field information exchange is carried out by converting the velocity into the form of absolute velocity. The two coordinate systems are as Figure 3 shown. The relative velocity in the rotating coordinate system is expressed as:

[0060]

[0061] In the formula, is the absolute velocity in the stationary coordinate system, is the rotational angular velocity of the wind turbine in the stationary coordinate system, and r is the spatial position vector.

[0062] For a wind turbine in the operating state, due to the large Reynolds number of the fluid (generally greater than 10^5), the flow is usually considered to be in a typical turbulent state. From a physical structure perspective, turbulence can be regarded as composed of vortices of different sizes superimposed on each other. Due to the interaction and random motion between vortices, turbulence has the characteristic of physical quantity pulsation. The unsteady mass equation and N - S (Navier - Stokes) equation can still be used to describe turbulence, and its air flow field model is as follows:

[0063]

[0064] Where: represents the air velocity; u, v, and w respectively represent the components of the air velocity in the x, y, and z directions, p represents the atmospheric pressure; υ represents the air dynamic viscosity.

[0065] For the convenience of numerical solution, the present invention uses the Reynolds - averaged equation to simplify the complex turbulent motion. According to the Reynolds - averaging method, the time - average value of any main control variable φ is defined as:

[0066]

[0067] In the formula, the superscript "-" represents the average over time. Using φ' to represent the pulsation value of the variable, the instantaneous value of the main control variable can be written as the sum of the average value and the pulsation value, expressed as:

[0068]

[0069] By introducing the idea of Reynolds averaging into the Navier - Stokes equation, the Reynolds - averaged N - S equation can be obtained; based on the incompressibility of the external flow field of the wind turbine, the viscosity coefficient of the fluid is a constant, and the energy equation is decoupled. In the Cartesian coordinate system, the mass equation of the three - dimensional air flow field is expressed as:

[0070]

[0071] In the formula, u i , u j (i, j = x, y, z) are the velocity components in the x, y, and z coordinate directions.

[0072] In this embodiment, determining the turbulent kinetic energy and the specific dissipation rate based on the air flow field model, and determining the derivatives of the turbulent kinetic energy and the specific dissipation rate with respect to time specifically include:

[0073] Constructing the mass equation of the three - dimensional air flow field in the Cartesian coordinate system:

[0074] where i, j = x, y, z; u' is the pulsating instantaneous value, ρ is the air density, and μ is the kinematic viscosity coefficient. represents the turbulent stress;

[0075] Due to the existence of the Reynolds stress term the number of unknowns is greater than the number of equations, resulting in the inability to close the system of equations. Therefore, a new turbulent equation needs to be introduced for solution. According to the eddy viscosity hypothesis proposed by Boussinesq, the turbulent viscosity coefficient μ t is introduced to relate the Reynolds stress and the mean velocity gradient as follows:

[0076] where: P t is the pressure caused by the pulsating velocity; where: k is the turbulent kinetic energy per unit mass of fluid;

[0077] δ ij is the Dirac function, and its expression is:

[0078]

[0079] For the external flow field of the wind turbine, it belongs to the typical high Reynolds number free flow condition. Near the blade surface, due to the viscous damping reducing the tangential velocity pulsation and the wall blocking the normal velocity, the turbulent flow development in the wind turbine blade surface area is not sufficient, so the Reynolds number is relatively low. The present invention uses the Shear Stress Transport (SST) k-ω turbulence model to close the Reynolds-averaged equations, solves the k-ω equation applicable to low Reynolds numbers in the near-wall region, solves the k-ε equation applicable to high Reynolds numbers in the outer boundary layer region and the free flow, and mixes these two turbulence models in the boundary layer.

[0080] μ t represents the turbulent viscosity coefficient. In the k-ω equation for low Reynolds numbers:

[0081]

[0082] In the k-ε equation for high Reynolds numbers:

[0083]

[0084] Determine the derivatives of the turbulent kinetic energy k and the specific dissipation rate ω with respect to time t:

[0085] where, P k represents the pressure caused by turbulence; ω is the specific dissipation rate; k is the turbulent kinetic energy; β * , σ k are closure constants; μ t is the turbulent viscosity coefficient.

[0086]

[0087] Wherein: F1 and F2 are auxiliary functions respectively, wherein:

[0088]

[0089] In this embodiment, the Lagrangian method is adopted to analyze the forces on the droplets in the gas-liquid two-phase flow field (this analysis process is prior art and will not be elaborated here), and the control equation for the water droplet phase in the air is constructed, specifically including:

[0090]

[0091] In the rectangular coordinate system, it is rewritten in matrix form:

[0092]

[0093] Wherein: K is the air-water droplet exchange coefficient, F x is the additional force caused by rotation;

[0094]

[0095] C D is the drag coefficient, R e represents the relative Reynolds number;

[0096]

[0097] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and not 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 the technical solutions of the present invention can be modified or equivalently replaced without departing from the spirit and scope of the technical solutions of the present invention, and they should all be covered within the scope of the claims of the present invention.

Claims

1. A method for determining the water droplet collision coefficient of a rotating fan blade based on the Lagrangian method, characterized in that: It includes the following steps: S1. Obtain the environmental parameters and the fan speed in the actual working condition, where the environmental parameters include the incoming flow wind speed, the droplet size in the air, and the atmospheric pressure; S2. Based on the environmental parameters and the fan speed in the actual working condition, simulate the rotation of the fan using the moving reference frame method and establish an air flow field model; S3. Determine the turbulent kinetic energy and the specific dissipation rate based on the air flow field model, and determine the derivatives of the turbulent kinetic energy and the specific dissipation rate with respect to time; S4. Analyze the forces acting on the droplets in the gas-liquid two-phase flow field using the Lagrangian method and construct the control equation for the water droplet phase in the air; S5. Input the derivatives of the turbulent kinetic energy and the specific dissipation rate with respect to time and the control equation for the water droplet phase into the DPM model in the FLUENT software to determine the distance between adjacent droplets at the release position and the distance between adjacent droplets hitting the fan blades; S6. Determine the droplet collision coefficient β: Where: dY represents the distance between adjacent water droplets at the release position, dL represents the distance between adjacent water droplets hitting the fan blade, β1 is the water droplet collision coefficient on the two-dimensional plane of the fan blade airfoil, and z1 and z2 are the distances of the positions where adjacent water droplets collide with the blade relative to the origin of the z-axis coordinate respectively.

2. The method for determining the water droplet collision coefficient of the rotating fan blade based on the Lagrange method according to claim 1, characterized in that: The air flow field model is: Wherein: represents the air velocity; u, v, and w respectively represent the components of the air velocity in the x, y, and z directions, p represents the atmospheric pressure; υ represents the air dynamic viscosity.

3. The method for determining the water droplet collision coefficient of the rotating fan blade based on the Lagrange method according to claim 2, characterized in that: Determining the turbulent kinetic energy and the specific dissipation rate based on the air flow field model and determining the derivatives of the turbulent kinetic energy and the specific dissipation rate with respect to time specifically include: Construct the mass equation of the three-dimensional air flow field in the Cartesian coordinate system: where i, j = x, y, z; ρ is the air density, μ is the kinematic viscosity coefficient, represents the turbulent stress; Where: P t is the pressure caused by the pulsating velocity; Where: k is the turbulent kinetic energy of the fluid per unit mass; δ ij is the Dirac function, and its expression is: μ t represents the turbulent viscosity coefficient, in the k-ω equation for low Reynolds numbers: In the high Reynolds number k-ε equation: Determine the derivatives of the turbulent kinetic energy k and the specific dissipation rate ω with respect to time t: Among them, P k represents the pressure caused by turbulence; ω is the specific dissipation rate; k is the turbulent kinetic energy; β * , σ k are closure constants; μ t is the turbulent viscosity coefficient. Where: F1 and F2 are auxiliary functions respectively, where:

4. The method for determining the water droplet collision coefficient of the rotating fan blade based on the Lagrange method according to claim 3, wherein: Constructing the control equation for the water droplet phase in the air specifically includes: In the rectangular coordinate system, rewrite it in matrix form: Where: K is the air-water droplet exchange coefficient, F x is the additional force caused by rotation; C D is the drag coefficient, R e represents the relative Reynolds number;

Images