A method for numerical simulation and distribution prediction of icing growth on a fan blade surface

Abstract

The present application relates to the field of fan blade deicing, and provides a fan blade surface icing growth numerical simulation and distribution prediction method, comprising: establishing a fan blade three-dimensional geometric model based on actual operation conditions and setting a rotating field flow field initial condition, constructing a rotating blade micro-domain cloud water flow field impact characteristic dynamic analysis based on an azimuth function, establishing a dynamic boundary slip model for blade surface liquid film formation and migration, coupling a thermodynamic phase change model and performing ice type dynamic growth calculation, performing iterative convergence judgment and icing form output based on geometric reconstruction, and performing aerodynamic and load analysis based on the final converged icing geometric model. The present application realizes dynamic reconstruction of the instantaneous impact characteristics of supercooled water droplets at different spatial positions of the rotating blade, can capture the non-uniform distribution of local liquid water content and water droplet diameter caused by terrain shielding or lifting effect when the blade rotates to different azimuth angles, and makes the input conditions of the icing simulation closer to the actual operation scene of the fan.

Description

Technical Field

[0001] This invention relates to the field of wind turbine blade de-icing technology, and in particular to a method for numerical simulation and distribution prediction of ice growth on the surface of wind turbine blades. Background Technology

[0002] Blade icing is a critical factor affecting the safe and efficient operation of wind turbines in cold regions. Especially in complex terrain conditions such as plateaus and mountains, the irregular growth of overflow ice can severely alter the aerodynamic characteristics of wind turbine blade airfoils, leading to increased power generation losses, blade load imbalances, and even turbine shutdowns. The formation and evolution of this problem are usually driven by the microphysical processes of supercooled water droplets impacting the blade surface. Among these factors, local micro-meteorological differences caused by blade rotation and the dynamic coupling between icing shape and thermodynamic parameters both affect the accuracy of icing pattern prediction.

[0003] In the field of wind turbine blade anti-icing technology, especially in predicting the icing morphology of large wind turbines in complex terrains such as plateaus and mountains, numerical simulation methods have become an important auxiliary tool for assessing icing risk and optimizing anti-icing strategies. Existing technologies typically rely on computational fluid dynamics and droplet impact models to predict icing growth on the blade surface by solving for the flow field distribution under uniform atmospheric environmental parameters. However, in practical applications, existing technologies lack corresponding coupled modeling mechanisms for the non-uniform distribution of localized supercooled water droplet concentration and size caused by terrain shading or lifting effects when the blade rotates to different azimuth angles, as well as the dynamic influence of ice shape changes on surface thermodynamic boundary conditions after icing growth. This leads to deviations between simulation results and the icing morphology under actual blade operating conditions. Summary of the Invention

[0004] The present invention aims to solve at least one of the problems existing in the prior art, and provides a method for numerical simulation and distribution prediction of ice growth on the surface of wind turbine blades.

[0005] This invention provides a method for numerical simulation and distribution prediction of icing growth on wind turbine blade surfaces, including: Step S1: Obtain the blade geometric parameters of the target wind turbine, the real-time operating parameters of the wind turbine unit, and the environmental meteorological parameters, construct a three-dimensional geometric model of the wind turbine blade, and set the initial flow field boundary conditions including the rotation domain. Step S2: Discretize the blades of the target wind turbine along the spanwise direction into multiple blade element micro-segments. Based on the real-time azimuth function of the blades, introduce the terrain disturbance correction factor to dynamically reconstruct the local subcooled water droplet impact characteristic parameters encountered by each blade element micro-segment at different azimuth angles during the rotation cycle. Step S3: Based on the impact characteristic parameters of local supercooled water droplets, a body-fitted coordinate system is constructed on the blade surface in the three-dimensional geometric model of the wind turbine blade, which grows with the ice shape. A dynamic boundary slip model considering the formation and migration of liquid film is established, taking into account centrifugal force, airflow shear stress and dynamic changes in ice surface roughness. The spatiotemporal distribution of liquid film thickness and liquid film flow velocity on the blade surface are calculated. Step S4: Based on the principle of thermodynamic phase change, divide the blade surface in the three-dimensional geometric model of the wind turbine blade into finite volume control units, solve the energy balance equation of each finite volume control unit, calculate the local icing rate, and determine the anisotropic direction of ice growth based on the vector sum of the liquid film flow direction and the local airflow direction, and generate the ice shape geometric data for the current time step. Step S5: Based on the ice-type geometry data, the computational mesh around the blade is updated using dynamic mesh technology. With the updated ice-type geometry as the boundary, steps S2 to S4 are repeated for iterative calculation until the ice mass change and the blade surface temperature field meet the convergence conditions. The three-dimensional blade ice morphology distribution data at the current moment is then output. Step S6: Perform aerodynamic and load analysis based on the finally converged icing geometry model to evaluate the impact of icing on wind turbine power and blade load.

[0006] Optionally, step S1 includes: Step S11: Obtain the actual airfoil data of the target wind turbine blades as blade geometric parameters, and establish a full-size three-dimensional geometric model including the target wind turbine blades, hub and near-wall boundary layer as the three-dimensional geometric model of the wind turbine blades; Step S12: Collect real-time operating parameters of the wind turbine, including rated speed and real-time pitch angle, as well as environmental meteorological parameters, including ambient temperature, ambient pressure, ambient wind speed, inflow liquid water content, and median volume diameter of water droplets. Step S13: In the computational fluid dynamics software, set up a rotating mechanical dynamic mesh region based on the three-dimensional geometric model of the wind turbine blade, define the blade surface in the three-dimensional geometric model of the wind turbine blade as a non-slip wall, set the rotational speed of the rotation domain, and consider the yaw error angle. Decompose the ambient wind speed vector into axial induced velocity and tangential induced velocity to obtain the initial flow field boundary conditions.

[0007] Optionally, step S2 includes: Step S21: Discretize the blades of the target wind turbine along the spanwise direction into multiple blade element segments. For each blade element segment and its current azimuth angle, calculate the relative composite velocity at the corresponding blade element segment based on the axial induction factor, tangential induction factor, and yaw error angle. Step S22: Based on the digital elevation model of the terrain where the target wind turbine is located, pre-calculate the changes in supercooled water droplet concentration caused by the terrain influence of the incoming flow at different azimuth angles, and generate a terrain disturbance correction factor; Step S23: Based on the relative synthesis velocity and terrain disturbance correction factor, combined with the incoming liquid water content, calculate the local liquid water content after considering the changes in terrain and relative velocity, and dynamically analyze the impact characteristics of the micro-domain cloud water flow field of the rotating blade. Step S24: In the normal direction of each control unit on the blade surface, a group of supercooled water droplets with representative diameters are launched. The diameter distribution of the supercooled water droplets is fitted with the Ross-Lammler distribution based on the measured median volume diameter of the water droplets. The initial velocity vector of the water droplets is set as the relative composite velocity vector. Step S25: Construct the control equation for water droplet motion based on the water droplet velocity to obtain the impact trajectory of the water droplet moving from the far field to the blade surface; solve the control equation for water droplet motion by numerical integration and count the amount of water droplet impact on each leaf element microsegment; Step S26: Determine the local collection coefficient based on the water droplet impact amount, local liquid water content, and relative synthesis velocity; calculate the local water droplet impact rate based on the local collection coefficient, local liquid water content, and relative synthesis velocity.

[0008] Optionally, the formula for calculating the relative synthesis rate is: ; in, For the first A leaf element microsegment at azimuth angle The relative synthesis rate at that location For ambient wind speed, The angle between the axis of the wind turbine nacelle and the direction of the incoming flow. The angular velocity of the wind turbine rotation. For the first The radius at a leaf element micro-segment As an axial inducing factor, Tangential induction factor; The formula for calculating the local liquid water content is: ; in, For the first A leaf element microsegment at azimuth angle The local liquid water content at that location The content of liquid water in the incoming flow. This is the compressibility correction factor. Azimuth The terrain disturbance correction factor at the location, and , Azimuth angle after taking into account the influence of terrain The actual local liquid water content of the incoming flow direction The original liquid water content of the environmental flow; The formula for calculating the local water droplet impact rate is: ; in, For the first A leaf element microsegment at azimuth angle The local water droplet impact rate at the location For the first A leaf element microsegment at azimuth angle The local collection coefficient at that location, and , The density of water, This represents the area of ​​the blade surface control unit.

[0009] Optionally, step S3 includes: Step S31: Based on the local water droplet impact rate, establish the governing equation for the liquid film flow to describe the relationship between the rate of change of liquid film mass over time and the water droplet impact rate, freezing rate, and evaporation rate. Step S32: Establish a liquid film momentum equation that considers the effects of centrifugal force, airflow shear stress and ice surface roughness, introduce a damping function based on rough Reynolds number, correct the near-wall viscous sublayer, and calculate the average flow velocity of the liquid film on the surface of the wind turbine blade. Step S33: Solve the governing equations and momentum equations of the liquid film flow simultaneously to obtain the spatiotemporal distribution data of the liquid film thickness on the blade surface.

[0010] Alternatively, the governing equation for liquid film flow can be expressed as: ; in, The density of water, The thickness of the liquid film on the surface of the wind turbine blades. Indicates time, For surface gradient operators, Let be the average flow velocity vector of the liquid film on the surface of the wind turbine blade. For local water droplet impact rate, For local icing rate, Evaporation / sublimation rate; The average velocity vector of the liquid film on the surface of the wind turbine blade The calculation formula is: ; in, For airflow shear stress, It is the centrifugal force vector. The area of ​​the blade surface control unit. The dynamic viscosity of water, For rough Reynolds number The damping function.

[0011] Optionally, step S4 includes: Step S41: Divide the blade surface in the three-dimensional geometric model of the wind turbine blade into finite volume control units, and establish an energy balance equation for each finite volume control unit, including water droplet impact heat, latent heat of phase change, convective heat transfer loss, evaporative heat dissipation, internal heat conduction of the blade, and aerodynamic heating. Step S42: Based on the energy balance equation, solve for the surface temperature of the wind turbine blades and calculate the local icing rate when the surface temperature of the wind turbine blades is below the freezing point. Step S43: Based on the liquid film flow direction and the local airflow direction, construct an anisotropic growth direction model, introduce the anisotropic growth coefficient determined by the local water film Froude number, and perform weighted synthesis of the surface normal and liquid film flow direction to determine the ice layer growth direction; Step S44: Based on the local icing rate and ice growth direction, calculate the displacement of the grid nodes on the wind turbine blade surface within the current time step, and generate the ice shape geometry data for the current time step.

[0012] Optionally, the formula for calculating the local icing rate when the surface temperature of the wind turbine blades is below the freezing point is: ; in, For local icing rate, The heat generated by the impact of water droplets pneumatic heating, For convective heat transfer losses, For evaporative heat dissipation, For heat conduction inside the fan blades This is the latent heat of freezing.

[0013] Alternatively, the direction of ice growth can be represented as follows: ; in, The direction of ice growth. Let be the unit vector normal to the surface of the wind turbine blade. Let be the unit vector in the direction of liquid film flow. The anisotropic growth coefficient and , For Froude number and , It is the acceleration due to gravity. For liquid film thickness, The average flow velocity of the liquid film on the surface of the wind turbine blades is denoted as .

[0014] Optionally, step S5 includes: Step S51: Based on the displacement of the mesh nodes on the surface of the wind turbine blade in the ice-type geometry data, update the computational mesh around the blade using the spring smoothing method or the local mesh reconstruction method to adapt the mesh to the updated ice-type geometry boundary. Step S52: Based on the updated ice-type geometric boundary, repeat steps S2 to S4 to calculate the updated flow field and icing rate; Step S53: Determine whether the change in icing mass within a single time step is less than a preset threshold and the average temperature field on the blade surface converges. If the conditions are met, stop the calculation; otherwise, repeat steps S2 to S4 to continue the iteration. Step S54: Output the three-dimensional blade icing morphology data at the current moment. The three-dimensional blade icing morphology data includes the distribution of ice thickness along the spanwise and tangential directions, ice density, and ice type characteristic parameters of key areas.

[0015] Compared with the prior art, the present invention has the following advantages: 1. By introducing an azimuth function and a terrain disturbance correction factor, the dynamic reconstruction of the instantaneous impact characteristics of supercooled water droplets at different spatial positions of the rotating blade is realized. This invention can capture the non-uniform distribution of local liquid water content and droplet diameter caused by terrain shading or lifting effect when the blade rotates to different azimuth angles, making the input conditions of icing simulation closer to the actual operating scenario of the wind turbine and improving the prediction accuracy of the icing initiation position and local growth rate. 2. In the calculation of liquid film migration on the blade surface, a slip boundary model considering centrifugal force, airflow shear stress and dynamic changes in ice surface roughness was established. This model can realistically reflect the asymmetric flow behavior of liquid film on the pressure and suction surfaces of the blade during icing. By introducing anisotropic growth coefficient to correct the direction of ice growth, the simulation problem of irregular growth of overflow ice towards the leeward side was effectively solved. 3. A geometric reconstruction technique with bidirectional iteration of ice shape growth and flow field calculation is adopted. The computational grid is regenerated and the flow field parameters are solved according to the updated ice shape at each time step. This allows the convective heat transfer coefficient and water droplet collection coefficient to be corrected in real time as the ice shape evolves, ensuring the thermodynamic balance and computational convergence in the long-term prediction of icing morphology. Attached Figure Description

[0016] One or more embodiments are illustrated by way of example with the corresponding pictures in the accompanying drawings. These illustrations do not constitute a limitation on the embodiments. Elements with the same reference numerals in the drawings are denoted as similar elements. Unless otherwise stated, the figures in the drawings are not to be limited by scale.

[0017] Figure 1 A flowchart illustrating a numerical simulation and distribution prediction method for ice growth on the surface of a wind turbine blade, provided as an embodiment of the present invention. Detailed Implementation

[0018] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the various embodiments of the present invention will be described in detail below with reference to the accompanying drawings. However, those skilled in the art will understand that many technical details are presented in the various embodiments of the present invention to facilitate a better understanding of the invention. However, the technical solutions claimed in the present invention can be implemented even without these technical details and with various variations and modifications based on the following embodiments. The division of the various embodiments below is for ease of description and should not constitute any limitation on the specific implementation of the present invention. The various embodiments can be combined with and referenced by each other without contradiction.

[0019] One embodiment of the present invention relates to a method for numerical simulation and distribution prediction of icing growth on the surface of wind turbine blades, the process of which is as follows: Figure 1 As shown, it includes steps S1 to S6.

[0020] Step S1: Obtain the blade geometric parameters of the target wind turbine, the real-time operating parameters of the wind turbine unit, and the environmental meteorological parameters, construct a three-dimensional geometric model of the wind turbine blade, and set the initial flow field boundary conditions including the rotation domain.

[0021] Specifically, step S1 is used to establish a three-dimensional geometric model of the wind turbine blades based on actual operating conditions and to set the initial conditions of the flow field in the rotating domain. For example, step S1 includes the following steps S11 to S13.

[0022] Step S11: Obtain the actual airfoil data of the target wind turbine blades as blade geometric parameters, and establish a full-size three-dimensional geometric model including the target wind turbine blades, hub and near-wall boundary layer as the three-dimensional geometric model of the wind turbine blades.

[0023] Specifically, step S11 is used to construct a full-size three-dimensional geometric model including the blade, hub and near-wall boundary layer based on the actual airfoil data of the target wind turbine blade, thereby obtaining the three-dimensional geometric model of the wind turbine blade.

[0024] Step S12: Collect real-time operating parameters of the wind turbine, including rated speed and real-time pitch angle, as well as environmental meteorological parameters, including ambient temperature, ambient pressure, ambient wind speed, inflow liquid water content, and median volume diameter of water droplets.

[0025] Specifically, step S12 is used to collect real-time operating parameters and environmental meteorological parameters of the wind turbine. The current real-time operating parameters of the wind turbine include: rated speed. Real-time pitch angle Environmental meteorological parameters include: ambient temperature. Environmental pressure Ambient wind speed Incoming liquid water content Median volume diameter of water droplets .

[0026] Step S13: In the computational fluid dynamics software, set up a rotating mechanical dynamic mesh region based on the three-dimensional geometric model of the wind turbine blade, define the blade surface in the three-dimensional geometric model of the wind turbine blade as a non-slip wall, set the rotational speed of the rotation domain, and consider the yaw error angle. Decompose the ambient wind speed vector into axial induced velocity and tangential induced velocity to obtain the initial flow field boundary conditions.

[0027] Specifically, in the computational fluid dynamics software, a moving mesh region for the rotating machinery is set, the blade surface is defined as a no-slip wall, and the rotational speed of the rotational domain is set to [value missing]. And taking into account the yaw error angle The environmental wind speed vector is decomposed into axial induced velocity and tangential induced velocity, laying the foundation for subsequent micro-domain flow field analysis.

[0028] Step S2: Discretize the blades of the target wind turbine along the spanwise direction into multiple blade element micro-segments. Based on the real-time azimuth function of the blades, introduce the terrain disturbance correction factor to dynamically reconstruct the local subcooled water droplet impact characteristic parameters encountered by each blade element micro-segment at different azimuth angles during the rotation cycle.

[0029] Specifically, step S2 is used to construct a dynamic analysis of the impact characteristics of the micro-domain cloud-water flow field of the rotating blade based on the azimuth angle function. For example, step S2 includes the following steps S21 to S26.

[0030] Step S21: Discretize the blades of the target wind turbine along the spanwise direction into multiple blade element segments. For each blade element segment and its current azimuth angle, calculate the relative composite velocity at the corresponding blade element segment based on the axial induction factor, tangential induction factor, and yaw error angle.

[0031] Specifically, the blade is discretized along the spanwise direction as follows: Each leaf element microsegment, targeting each leaf element microsegment and its current azimuth angle It is necessary to calculate the relative synthesis rate at this leaf nutrient segment, where, The time during one rotation of the blade The angle of change. The relative composite velocity is a fundamental parameter determining the impact kinetic energy of water droplets and the local collection coefficient; its calculation formula comprehensively considers the axial induction factor. Tangential inducing factor The changes in velocity components caused by terrain provide accurate velocity boundary conditions for subsequent water droplet impact calculations. For example, the formula for calculating the relative composite velocity is: ; in, For the first A leaf element microsegment at azimuth angle The relative synthesis rate at that location For ambient wind speed, The angle between the axis of the wind turbine nacelle and the direction of the incoming flow. The angular velocity of the wind turbine rotation. For the first The radius at a leaf element micro-segment As an axial inducing factor, It is a tangential inducing factor.

[0032] Step S22: Based on the digital elevation model of the terrain where the target wind turbine is located, pre-calculate the changes in supercooled water droplet concentration caused by the terrain influence of the incoming flow at different azimuth angles, and generate a terrain disturbance correction factor.

[0033] Specifically, since wind turbines are typically located in high-altitude and mountainous areas, the undulating terrain can lead to uneven distribution of subcooled water droplet concentration. Therefore, this embodiment introduces a terrain disturbance correction factor. This terrain disturbance correction factor is based on a digital elevation model that pre-calculates the impact of terrain lifting or shading on incoming flows at different azimuth angles. For example, the formula for calculating the terrain disturbance correction factor is: ; in, Azimuth Terrain disturbance correction factor at the location. Azimuth angle after taking into account the influence of terrain The actual local liquid water content of the directional flow; The original liquid water content of the environment.

[0034] Step S23: Based on the relative synthesis velocity and terrain disturbance correction factor, combined with the incoming liquid water content, calculate the local liquid water content after considering the changes in terrain and relative velocity, and dynamically analyze the impact characteristics of the micro-domain cloud water flow field of the rotating blade.

[0035] Specifically, taking into account topography and changes in relative synthesis rate, and combining leaf element microfractions In azimuth Calculate the local liquid water content based on the actual concentration of supercooled water droplets encountered. This method corrects errors introduced by the traditional assumption of uniform flow, making water droplet impact calculations more consistent with the realities of mountain wind farms. For example, the formula for calculating localized liquid water content is: ; In the formula, For the first A leaf element microsegment at azimuth angle The local liquid water content at that location The liquid water content of the incoming flow is the liquid water content of the environment. This is the compressibility correction factor. In particular, the value of k can be determined based on the physical magnitude. For example, under normal wind speed and blade tip Mach number below 0.3, the compressibility effect can be ignored, and k=0 can be directly taken; while when the compressibility effect cannot be ignored, the value of k is in the range of 0.1-0.3, that is, starting from 0.1 and not exceeding 0.3.

[0036] The droplet impact rate refers to the mass of supercooled water droplets impacting a unit area of ​​the blade surface per unit time. It is used as input for the mass source term in the subsequent liquid film continuity equation. The droplet impact rate cannot be directly obtained through measurement; it requires information based on the local collection coefficient of the water droplets. The local liquid water content and relative synthesis rate are jointly determined. Specifically, this can be achieved by using an improved droplet trajectory tracking algorithm based on the Lagrange method, tailored to each micro-segment of the leaf. and azimuth Perform the calculation.

[0037] Step S24: In the normal direction of each control unit on the blade surface, launch a set of representative diameters. The diameter distribution of supercooled water droplets is based on the measured median volume diameter of the droplets. The Rosin-Rammler distribution was used for fitting, and the initial velocity vector of the water droplet was set as the relative composite velocity vector.

[0038] Specifically, the control unit on the blade surface refers to the finite volume control unit obtained by dividing the blade surface in the three-dimensional geometric model of the wind turbine blade. The initial velocity vector of the water droplet. Let it be the relative resultant velocity vector Right now That is, there is no velocity slip between the water droplet and the airflow (assuming that the water droplet follows the airflow until it approaches the wall).

[0039] Step S25: Construct the control equation for water droplet motion based on the water droplet velocity to obtain the impact trajectory of the water droplet moving from the far field to the blade surface; solve the control equation for water droplet motion by numerical integration and count the amount of water droplet impact on each leaf element microsegment.

[0040] Specifically, water droplets in an airflow are primarily affected by drag and gravity, neglecting secondary factors such as lift and pressure gradient force. The acceleration of the water droplet is equal to the vector sum of the acceleration caused by the airflow drag and the acceleration due to gravity, with the drag coefficient being... The following empirical formula can be used to calculate it: ,in, Let be the relative Reynolds number of the water droplet.

[0041] The governing equations for the motion of a water droplet are used to calculate the trajectory of the droplet as it moves from the far field to the surface of the blade. The governing equations for the motion of a water droplet can be expressed as follows: ; in, Let be the velocity vector of the water droplet. Aerodynamic viscosity, The density of water; The local airflow velocity vector, the specific value of which can be determined by the flow field containing the rotating domain in step S1; It is the gravitational acceleration vector. The diameter of the water droplet. The relative Reynolds number of the water droplet; The above control equations for the motion of water droplets are solved by numerical integration to obtain the position and angle of the water droplets impacting the blade surface; based on the distribution of impact points, the amount of water droplet impact on each blade surface control unit is statistically analyzed.

[0042] Step S26: Determine the local collection coefficient based on the water droplet impact amount, local liquid water content, and relative synthesis velocity; calculate the local water droplet impact rate based on the local collection coefficient, local liquid water content, and relative synthesis velocity.

[0043] Specifically, define the local collection coefficient. This is the ratio of the actual water droplet mass flux received by the control unit corresponding to the leaf element microsegment to the incoming water droplet mass flux based on the relative synthesis velocity. The local collection coefficient represents the efficiency of water droplets impacting the wall due to inertia deviating from the streamline. The larger the value, the easier it is for water droplets to be collected at that location. The local collection coefficient can be used to convert flow field information into icing mass input. For example, the formula for calculating the local collection coefficient is: ; Here, the control unit refers to the first The control unit corresponding to each leaf element microsegment The area of ​​the blade surface control unit. This is the density of water.

[0044] Local water droplet impact rate The mass source term is used as the governing equation for the liquid film flow in step S3. An example is the local droplet impact rate. The calculation formula is: ; in, For the first A leaf element microsegment at azimuth angle The local water droplet impact rate at a given location.

[0045] The above calculations can be used to obtain different azimuth angles. Different leaf element segments Real-time water droplet impact rate for each surface control unit.

[0046] Step S3: Based on the impact characteristic parameters of local supercooled water droplets, construct a body-fitted coordinate system on the blade surface in the three-dimensional geometric model of the wind turbine blade, which grows with the ice shape. Establish a dynamic boundary slip model for liquid film formation and migration that considers centrifugal force, airflow shear stress, and dynamic changes in ice surface roughness. Calculate the spatiotemporal distribution of liquid film thickness and liquid film flow velocity on the blade surface.

[0047] Specifically, step S3 is used to establish a dynamic boundary slip model for the formation and migration of liquid film on the blade surface. For example, step S3 includes steps S31 to S33.

[0048] Step S31: Based on the local water droplet impact rate, establish the governing equation for the liquid film flow to describe the relationship between the rate of change of liquid film mass over time and the water droplet impact rate, freezing rate, and evaporation rate.

[0049] Specifically, a body-fitted coordinate system that grows with the ice type is constructed on the blade surface. The local water droplet impact rate is calculated based on step S2. Therefore, a governing equation for liquid film flow needs to be established. The rate of change of liquid film mass over time is equal to the local droplet impact rate minus the local freezing rate and evaporation rate. This equation is used to quantitatively track the accumulation and consumption of liquid water on the blade surface, providing the water film distribution for subsequent freezing calculations. For example, the governing equation for liquid film flow is expressed as: ; in, The density of water, The thickness of the liquid film on the surface of the wind turbine blades. Indicates time, For surface gradient operators, This is the average velocity vector of the liquid film, which is the average flow velocity vector of the liquid film on the surface of the wind turbine blades. For local water droplet impact rate, For local icing rate, This represents the evaporation / sublimation rate.

[0050] Step S32: Establish a liquid film momentum equation that considers the effects of centrifugal force, airflow shear stress, and ice surface roughness. Introduce a damping function based on the rough Reynolds number to correct the near-wall viscous sublayer and calculate the average flow velocity of the liquid film on the wind turbine blade surface.

[0051] Specifically, the average velocity of the liquid film, i.e., the average flow velocity of the liquid film on the wind turbine blade surface, is jointly determined by the airflow shear stress, centrifugal force, and viscous subsurface damping. The damping function is dynamically adjusted according to the surface roughness of the ice layer, thus more realistically simulating the flow of the water film on a rough ice surface. This allows for accurate calculation of the non-uniform migration of the water film on the blade surface (pressure and suction sides), directly affecting the morphology of overflow ice. The average flow velocity vector of the liquid film on the wind turbine blade surface... The calculation formula is the analytical solution of the liquid film momentum equation under reasonable simplification. For example, the average flow velocity vector of the liquid film on the surface of the wind turbine blade... The calculation formula is: ; in, For airflow shear stress, It is the centrifugal force vector. The area of ​​the blade surface control unit. The dynamic viscosity of water, The liquid film thickness refers to the thickness of the liquid film on the surface of the wind turbine blades. For rough Reynolds number The damping function.

[0052] Step S33: Solve the governing equations and momentum equations of the liquid film flow simultaneously to obtain the spatiotemporal distribution data of the liquid film thickness on the blade surface.

[0053] Step S4: Based on the principle of thermodynamic phase change, divide the blade surface in the three-dimensional geometric model of the wind turbine blade into finite volume control units, solve the energy balance equation of each finite volume control unit, calculate the local icing rate, and determine the anisotropic direction of ice growth based on the vector sum of the liquid film flow direction and the local airflow direction, and generate the ice shape geometric data for the current time step.

[0054] Specifically, step S4 is used to couple the thermodynamic phase transition model with the dynamic growth calculation of ice. For example, step S4 includes steps S41 to S44.

[0055] Step S41: Divide the blade surface in the three-dimensional geometric model of the wind turbine blade into finite volume control units, and establish an energy balance equation for each finite volume control unit, including water droplet impact heat, latent heat of phase change, convective heat transfer loss, evaporative heat dissipation, internal heat conduction of the blade, and aerodynamic heating.

[0056] Specifically, the icing process is essentially a thermodynamic phase transition process, and the energy balance determines the ice growth rate and ice type. The energy terms include: heat flux from water droplet impact, latent heat released during phase transition, convective heat transfer losses, evaporative heat dissipation losses, internal thermal conductivity of the blade composite material, and aerodynamic heating. The energy balance equation can be expressed as: Heat flux from water droplet impact + Aerodynamic heating heat = Convective heat transfer losses + Evaporative heat dissipation losses + Internal thermal conductivity losses of the blade composite material + Latent heat released during icing.

[0057] Step S42: Based on the energy balance equation, solve for the surface temperature of the wind turbine blades and calculate the local icing rate when the surface temperature of the wind turbine blades is below the freezing point.

[0058] Specifically, the energy term in the energy balance equation implicitly contains the surface temperature of the wind turbine blades, and the surface temperature of the wind turbine blades can be obtained by solving the energy balance equation and related existing technologies. Introducing an interface temperature slip criterion: If the surface temperature of the fan blades... Below freezing point If the temperature is below the freezing point, icing is considered to have occurred. The amount of icing is determined by excess cooling. For example, the formula for calculating the local icing rate when the surface temperature of a fan blade is below the freezing point is: ; in, For local icing rate, Impact heat refers to the heat flux generated by the impact of water droplets. Pneumatic heating refers to the heat generated by pneumatic heating. For convective heat transfer losses, This refers to evaporative heat loss, or evaporative heat dissipation loss. This refers to the internal heat conduction of the wind turbine blades, specifically the heat loss within the composite material of the blades. The latent heat of freezing is the latent heat released when ice forms.

[0059] Step S43: Based on the liquid film flow direction and the local airflow direction, construct an anisotropic growth direction model, introduce the anisotropic growth coefficient determined by the local water film Froude number, and perform weighted synthesis of the surface normal and liquid film flow direction to determine the ice layer growth direction.

[0060] Specifically, the calculated local icing rate This will lead to changes in the blade surface geometry. By employing a bidirectional iterative geometry update technique, the ice growth direction is no longer simply along the normal, but is determined by the vector sum of the liquid film flow direction and the local airflow direction, to simulate the physical phenomenon of overflow ice growing towards the leeward side. For example, the ice growth direction is represented as: ; in, The direction of ice growth. Let be the unit vector normal to the surface of the wind turbine blade. The unit vector represents the direction of liquid film flow. The anisotropic growth coefficient is 0 ≤ ≤1, By Froude Decision and , Among them, Froude number A dimensionless number characterizing the ratio of inertial force to gravity in liquid film flow, used to determine the state of liquid film flow; The average flow velocity of the liquid film on the blade surface is determined by the combined effects of airflow shear stress, centrifugal force, and viscous damping. The value of is the average flow velocity vector of the liquid film on the surface of the wind turbine blades. Size; It is the acceleration due to gravity; The liquid film thickness refers to the thickness of the liquid film on the surface of the wind turbine blade, and it is also the thickness of the water film covering the blade surface in the normal direction.

[0061] Step S44: Based on the local icing rate and ice growth direction, calculate the displacement of the grid nodes on the wind turbine blade surface within the current time step, and generate the ice shape geometry data for the current time step.

[0062] Specifically, the ice geometry data for the current time step includes the calculated displacement of the grid nodes on the wind turbine blade surface within the current time step. The magnitude of the displacement of the grid nodes on the wind turbine blade surface within the current time step can be obtained by calculating the product of the local icing rate and the current time step, and the direction of displacement is the direction of ice growth.

[0063] Step S5: Based on the ice-type geometry data, the computational mesh around the blade is updated using dynamic meshing technology. With the updated ice-type geometry as the boundary, steps S2 to S4 are repeated for iterative calculation until the change in icing mass and the blade surface temperature field meet the convergence conditions. The three-dimensional blade icing morphology distribution data at the current moment is then output.

[0064] Specifically, step S5 is used to perform iterative convergence judgment and output the icing morphology based on geometric reconstruction. For example, step S5 includes the following steps S51 to S54.

[0065] Step S51: Based on the displacement of the grid nodes on the surface of the wind turbine blade in the ice-type geometry data, update the computational grid around the blade using the spring smoothing method or the local grid reconstruction method to adapt the grid to the updated ice-type geometry boundary.

[0066] Specifically, due to the change in ice shape, the original flow field mesh for clean blades is no longer applicable and must be regenerated. Therefore, in step S51, the computational mesh around the blades can be updated using spring smoothing or local reconstruction methods based on the displacement of the mesh nodes on the wind turbine blade surface calculated in step S4, thereby obtaining the updated ice shape geometric boundary, i.e., the new ice shape geometric boundary.

[0067] Step S52: Based on the updated ice-type geometric boundary, repeat steps S2 to S4 to calculate the updated flow field and icing rate.

[0068] Step S53: Determine whether the change in icing mass within a single time step is less than a preset threshold and the average temperature field on the blade surface converges. If the conditions are met, stop the calculation; otherwise, repeat steps S2 to S4 to continue the iteration.

[0069] Step S54: Output the three-dimensional blade icing morphology data at the current moment. The three-dimensional blade icing morphology data includes the distribution of ice thickness along the spanwise and tangential directions, ice density, and ice type characteristic parameters of key areas.

[0070] Specifically, using the new ice-shaped geometry as the boundary, steps S2 to S4 are repeated until a single time step is reached. The change in icing mass within the area is less than a set threshold. Furthermore, the average temperature field on the blade surface converges, outputting the current moment. The three-dimensional blade icing morphology data includes the distribution of ice thickness along the spanwise and tangential directions, ice density, and ice type characteristic parameters of key areas.

[0071] Step S6: Perform aerodynamic and load analysis based on the finally converged icing geometry model to evaluate the impact of icing on wind turbine power and blade load.

[0072] Specifically, step S6 can construct an aerodynamic and load database after icing, and use the finally converged icing geometric model, i.e. the final obtained three-dimensional geometric model of the iced wind turbine blade, to perform aeroelastic calculations based on existing aeroelastic simulation methods, calculate the aerodynamic performance attenuation coefficient of the blade after icing, and then combine it with relevant wind turbine control system strategies to predict the power loss and unbalanced load caused by icing, providing data support for subsequent de-icing start-up control or wind turbine shutdown protection.

[0073] The working principle of the numerical simulation and distribution prediction method for wind turbine blade surface icing growth provided in this invention is as follows: By establishing a transient model of the local micro-meteorological field during blade rotation, and combining it with a bidirectional iterative algorithm for surface liquid film migration and ice layer growth, real-time feedback of geometric changes on flow field and thermodynamic parameters is introduced into the icing morphology prediction, thereby enhancing the simulation accuracy of complex overflow ice morphology; the spectrum and impact velocity of supercooled water droplets encountered by each micro-segment of the blade are dynamically analyzed using an azimuth function, while a terrain disturbance correction factor is introduced to correct the local liquid water content, forcing the icing growth process to be integrated with surface water film flow, convective heat transfer, and internal blade processes. The heat conduction of the blade reaches energy and mass balance within each time step, thereby improving the prediction accuracy of irregular icing distribution at the leading edge and maximum thickness of the wind turbine blade under asymmetric meteorological input scenarios caused by blade rotation. By discretizing the blade spanwise into several blade element micro-segments and coupling rotational dynamics, a water droplet impact model based on local synthesis velocity is established, and an interfacial thermodynamic slip boundary algorithm is introduced, enabling the model to simultaneously capture the driving effect of centrifugal force on liquid film migration and the inhibitory effect of ice roughness on aerodynamic heat exchange. This avoids the distortion of icing rate calculation under traditional static assumptions while preserving the anisotropic details of ice growth.

[0074] In the numerical simulation and distribution prediction method for icing growth on wind turbine blade surfaces provided by this invention, local meteorological parameters on the blade element micro-segments are reconstructed based on real-time azimuth angles and topographic digital elevation models. The synthesized velocity is vector-synthesized by considering axial and tangential induction factors and yaw error angles, while the liquid water content is corrected based on the obstruction or lifting effect of topography on clouds and fog. These dynamic parameters are used to solve the Lagrange droplet trajectory or the Eulerian multiphase flow model, thereby quantifying the droplet collection coefficient and impact energy on each blade element micro-segment. The core of icing growth prediction lies in upgrading the traditional unidirectional flow-solid-thermal decoupling process to a bidirectional iterative coupling: firstly, within each time step, the continuity and momentum equations of the surface liquid film are solved based on the current ice-type geometry; the liquid film velocity is solved using an improved Van... The Driest damping function reflects the impediment of near-wall flow by ice roughness; the icing rate is then calculated using the energy balance equation, where the balance between latent heat of icing and various heat fluxes is automatically controlled by the interface temperature slip criterion; finally, the ice growth direction is determined based on the vector sum of the liquid film flow direction and the local airflow direction, and the ice geometry is updated using dynamic mesh technology. The above steps are then repeated with the new ice geometry as the boundary until convergence. The core of the bidirectional iterative mechanism is to update the liquid film velocity through surface roughness and control the ice extension direction through the anisotropic growth coefficient, while penalizing the error accumulation caused by the failure of ice shape changes to be fed back to the flow field in traditional methods; each update of the ice geometry forces the recalculation of the local convective heat transfer coefficient and the droplet collection coefficient: this feedback can effectively correlate the interaction between the ice shape and the local microclimate without the need for additional experimental calibration.

[0075] Through the aforementioned dynamic coupling and geometric reconstruction mechanism, the numerical simulation and distribution prediction method for ice growth on wind turbine blades provided in this invention can effectively solve the problem of simulating irregular overflow ice growth on large wind turbine blades in complex terrains of plateaus and mountains. This avoids significant errors caused by traditional uniform meteorological input models due to wind speed shear, uneven cloud and fog distribution, and fluctuations in supercooled water droplet concentration experienced by the blades during rotation. The dynamic coupling mechanism, by establishing an azimuth-dependent local micro-meteorological field and a bidirectional iteration of ice type-flow field, can restore physically meaningful structural patterns during icing. At the same time, it avoids blurring local thermodynamic details under fixed geometric assumptions, making the prediction results more adaptable and reliable to meteorological disturbances in actual operating conditions. The preprocessed three-dimensional blade geometric model, real-time meteorological parameters, and optimized iterative algorithm are ultimately used to output the distribution of ice thickness along the spanwise and chordwise directions, ice density, and key ice type characteristic parameters, thereby providing a reliable basis for the design of wind turbine anti-icing and de-icing systems, power loss assessment, and operation strategy adjustment. The entire process can be achieved through automated mesh reconstruction and convergence determination, ensuring that the prediction results retain ice type details while improving practicality in complex mountain wind farm scenarios.

[0076] Those skilled in the art will understand that the above embodiments are specific implementations of the present invention, and in practical applications, various changes can be made in form and detail without departing from the spirit and scope of the present invention.

Claims

1. A numerical simulation and distribution prediction method for ice growth on the surface of wind turbine blades, characterized in that, include: Step S1: Obtain the blade geometric parameters of the target wind turbine, the real-time operating parameters of the wind turbine unit, and the environmental meteorological parameters, construct a three-dimensional geometric model of the wind turbine blade, and set the initial flow field boundary conditions including the rotation domain. Step S2: Discretize the blades of the target wind turbine along the spanwise direction into multiple blade element micro-segments. Based on the real-time azimuth function of the blades, introduce the terrain disturbance correction factor to dynamically reconstruct the local subcooled water droplet impact characteristic parameters encountered by each blade element micro-segment at different azimuth angles during the rotation cycle. Step S3: Based on the impact characteristic parameters of local supercooled water droplets, a body-fitted coordinate system is constructed on the blade surface in the three-dimensional geometric model of the wind turbine blade, which grows with the ice shape. A dynamic boundary slip model considering the formation and migration of liquid film is established, taking into account centrifugal force, airflow shear stress and dynamic changes in ice surface roughness. The spatiotemporal distribution of liquid film thickness and liquid film flow velocity on the blade surface are calculated. Step S4: Based on the principle of thermodynamic phase change, divide the blade surface in the three-dimensional geometric model of the wind turbine blade into finite volume control units, solve the energy balance equation of each finite volume control unit, calculate the local icing rate, and determine the anisotropic direction of ice growth based on the vector sum of the liquid film flow direction and the local airflow direction, and generate the ice shape geometric data for the current time step. Step S5: Based on the ice-type geometry data, the computational mesh around the blade is updated using dynamic mesh technology. With the updated ice-type geometry as the boundary, steps S2 to S4 are repeated for iterative calculation until the ice mass change and the blade surface temperature field meet the convergence conditions. The three-dimensional blade ice morphology distribution data at the current moment is then output. Step S6: Perform aerodynamic and load analysis based on the finally converged icing geometry model to evaluate the impact of icing on wind turbine power and blade load.

2. The method for numerical simulation and distribution prediction of ice growth on wind turbine blade surface according to claim 1, characterized in that, Step S1 includes: Step S11: Obtain the actual airfoil data of the target wind turbine blades as blade geometric parameters, and establish a full-size three-dimensional geometric model including the target wind turbine blades, hub and near-wall boundary layer as the three-dimensional geometric model of the wind turbine blades; Step S12: Collect real-time operating parameters of the wind turbine, including rated speed and real-time pitch angle, as well as environmental meteorological parameters, including ambient temperature, ambient pressure, ambient wind speed, inflow liquid water content, and median volume diameter of water droplets. Step S13: In the computational fluid dynamics software, set up a rotating mechanical dynamic mesh region based on the three-dimensional geometric model of the wind turbine blade, define the blade surface in the three-dimensional geometric model of the wind turbine blade as a non-slip wall, set the rotational speed of the rotation domain, and consider the yaw error angle. Decompose the ambient wind speed vector into axial induced velocity and tangential induced velocity to obtain the initial flow field boundary conditions.

3. The method for numerical simulation and distribution prediction of ice growth on wind turbine blade surface according to claim 1, characterized in that, Step S2 includes: Step S21: Discretize the blades of the target wind turbine along the spanwise direction into multiple blade element segments. For each blade element segment and its current azimuth angle, calculate the relative composite velocity at the corresponding blade element segment based on the axial induction factor, tangential induction factor, and yaw error angle. Step S22: Based on the digital elevation model of the terrain where the target wind turbine is located, pre-calculate the changes in supercooled water droplet concentration caused by the terrain influence of the incoming flow at different azimuth angles, and generate a terrain disturbance correction factor; Step S23: Based on the relative synthesis velocity and terrain disturbance correction factor, combined with the incoming liquid water content, calculate the local liquid water content after considering the changes in terrain and relative velocity, and dynamically analyze the impact characteristics of the micro-domain cloud water flow field of the rotating blade. Step S24: In the normal direction of each control unit on the blade surface, a group of supercooled water droplets with representative diameters are launched. The diameter distribution of the supercooled water droplets is fitted with the Ross-Lammler distribution based on the measured median volume diameter of the water droplets. The initial velocity vector of the water droplets is set as the relative composite velocity vector. Step S25: Construct the control equation for water droplet motion based on the water droplet velocity to obtain the impact trajectory of the water droplet moving from the far field to the blade surface; solve the control equation for water droplet motion by numerical integration and count the amount of water droplet impact on each leaf element microsegment; Step S26: Determine the local collection coefficient based on the water droplet impact amount, local liquid water content, and relative synthesis velocity; calculate the local water droplet impact rate based on the local collection coefficient, local liquid water content, and relative synthesis velocity.

4. The method for numerical simulation and distribution prediction of ice growth on wind turbine blade surface according to claim 3, characterized in that, The formula for calculating the relative synthesis rate is: ； in, For the first A leaf element microsegment at azimuth angle The relative synthesis rate at that location For ambient wind speed, The angle between the axis of the wind turbine nacelle and the direction of the incoming flow. The angular velocity of the wind turbine rotation. For the first The radius at a leaf element micro-segment As an axial inducing factor, Tangential induction factor; The formula for calculating the local liquid water content is: ； in, For the first A leaf element microsegment at azimuth angle The local liquid water content at that location The content of liquid water in the incoming flow. This is the compressibility correction factor. Azimuth The terrain disturbance correction factor at the location, and , Azimuth angle after taking into account the influence of terrain The actual local liquid water content of the incoming flow direction The original liquid water content of the environmental flow; The formula for calculating the local water droplet impact rate is: ； in, For the first A leaf element microsegment at azimuth angle The local water droplet impact rate at the location For the first A leaf element microsegment at azimuth angle The local collection coefficient at that location, and , The density of water, This represents the area of ​​the blade surface control unit.

5. The method for numerical simulation and distribution prediction of ice growth on wind turbine blade surface according to claim 1, characterized in that, Step S3 includes: Step S31: Based on the local water droplet impact rate, establish the governing equation for the liquid film flow to describe the relationship between the rate of change of liquid film mass over time and the water droplet impact rate, freezing rate, and evaporation rate. Step S32: Establish a liquid film momentum equation that considers the effects of centrifugal force, airflow shear stress and ice surface roughness, introduce a damping function based on rough Reynolds number, correct the near-wall viscous sublayer, and calculate the average flow velocity of the liquid film on the surface of the wind turbine blade. Step S33: Solve the governing equations and momentum equations of the liquid film flow simultaneously to obtain the spatiotemporal distribution data of the liquid film thickness on the blade surface.

6. The method for numerical simulation and distribution prediction of ice growth on wind turbine blade surface according to claim 5, characterized in that, The governing equation for liquid film flow is expressed as: ； in, The density of water, The thickness of the liquid film on the surface of the wind turbine blades. Indicates time, For surface gradient operators, Let be the average flow velocity vector of the liquid film on the surface of the wind turbine blade. For local water droplet impact rate, For local icing rate, Evaporation / sublimation rate; The average velocity vector of the liquid film on the surface of the wind turbine blade The calculation formula is: ； in, For airflow shear stress, It is the centrifugal force vector. The area of ​​the blade surface control unit. The dynamic viscosity of water, For rough Reynolds number The damping function.

7. The method for numerical simulation and distribution prediction of ice growth on wind turbine blade surface according to claim 1, characterized in that, Step S4 includes: Step S41: Divide the blade surface in the three-dimensional geometric model of the wind turbine blade into finite volume control units, and establish an energy balance equation for each finite volume control unit, including water droplet impact heat, latent heat of phase change, convective heat transfer loss, evaporative heat dissipation, internal heat conduction of the blade, and aerodynamic heating. Step S42: Based on the energy balance equation, solve for the surface temperature of the wind turbine blades and calculate the local icing rate when the surface temperature of the wind turbine blades is below the freezing point. Step S43: Based on the liquid film flow direction and the local airflow direction, construct an anisotropic growth direction model, introduce the anisotropic growth coefficient determined by the local water film Froude number, and perform weighted synthesis of the surface normal and liquid film flow direction to determine the ice layer growth direction; Step S44: Based on the local icing rate and ice growth direction, calculate the displacement of the grid nodes on the wind turbine blade surface within the current time step, and generate the ice shape geometry data for the current time step.

8. The method for numerical simulation and distribution prediction of ice growth on wind turbine blade surface according to claim 7, characterized in that, The formula for calculating the local icing rate when the surface temperature of a wind turbine blade is below the freezing point is: ； in, For local icing rate, The heat generated by the impact of water droplets pneumatic heating, For convective heat transfer losses, For evaporative heat dissipation, For heat conduction inside the fan blades This is the latent heat of freezing.

9. The method for numerical simulation and distribution prediction of ice growth on wind turbine blade surface according to claim 7, characterized in that, The direction of ice growth is represented as follows: ； in, The direction of ice growth. Let be the unit vector normal to the surface of the wind turbine blade. Let be the unit vector in the direction of liquid film flow. The anisotropic growth coefficient and , For Froude number and , It is the acceleration due to gravity. For liquid film thickness, The average flow velocity of the liquid film on the surface of the wind turbine blades is denoted as .

10. The method for numerical simulation and distribution prediction of ice growth on wind turbine blade surface according to claim 1, characterized in that, Step S5 includes: Step S51: Based on the displacement of the mesh nodes on the surface of the wind turbine blade in the ice-type geometry data, update the computational mesh around the blade using the spring smoothing method or the local mesh reconstruction method to adapt the mesh to the updated ice-type geometry boundary. Step S52: Based on the updated ice-type geometric boundary, repeat steps S2 to S4 to calculate the updated flow field and icing rate; Step S53: Determine whether the change in icing mass within a single time step is less than a preset threshold and the average temperature field on the blade surface converges. If the conditions are met, stop the calculation; otherwise, repeat steps S2 to S4 to continue the iteration. Step S54: Output the three-dimensional blade icing morphology data at the current moment. The three-dimensional blade icing morphology data includes the distribution of ice thickness along the spanwise and tangential directions, ice density, and ice type characteristic parameters of key areas.

Images