A method and system for solving the coupling of water film flow-heat transfer-phase change in the icing process of a wind turbine blade

Abstract

The application discloses a wind turbine blade icing process water film flow-heat-phase change coupling solving method and system, and belongs to the technical field of wind power generation. The method comprises the following steps: by constructing a water film flow-heat-phase change coupling control equation group under a three-dimensional curved surface coordinate system, fully considering the influence of centrifugal force and Coriolis force caused by blade rotation, realizing strong coupling solving of water film flow and phase change process; adopting a partition adaptive time step strategy, combining a dynamic grid updating technology, and realizing real-time linkage of water film flow, heat transfer and phase change through a frozen proportion coefficient feedback mechanism. The application solves the problem that precision and efficiency are difficult to be considered in the prior art, improves the authenticity of ice shape simulation, and provides a precise and efficient numerical simulation scheme for wind turbine blade icing protection and related research.

Description

Technical Field

[0001] This invention relates to the fields of wind power generation technology and computational fluid dynamics, specifically to a method and system for solving the coupled water film flow-heat transfer-phase change process during the icing of wind turbine blades. Background Technology

[0002] When wind turbines operate in cold and humid climates, the blade surface is prone to icing. The icing process is essentially a complex physical process in which supercooled water droplets impact the blades, forming a water film. During the flow of this water film, heat transfer and phase change (freezing) occur simultaneously. Accurate simulation of this process is of great significance for predicting ice shape distribution, assessing the impact of icing, and designing anti-icing and de-icing systems.

[0003] In the prior art, there are already disclosed methods for icing simulation, such as the invention with publication number CN120724618B, entitled "A Machine Learning-Based Method for Predicting Icing on Aero-engine Rotating Blades." This invention discloses an icing prediction and simulation method that calculates the icing process by constructing a water film control equation and combining it with thermodynamic equilibrium, and can reflect the icing-related physical processes to a certain extent. However, even based on this type of prior art, current numerical simulations of icing on wind turbine blades still have the following significant shortcomings:

[0004] 1. The simplified description of water film flow characteristics makes it difficult to reflect the actual icing morphology: Most existing models still simplify the treatment of water film and do not fully consider the flow, overflow and redistribution of water film on the three-dimensional complex curved surface of wind turbine blades. It is difficult to accurately depict typical ice shapes such as ice corners and double-cornered ice formed by water film migration. The simulation results deviate significantly from the actual icing morphology.

[0005] 2. Imperfect coupling mechanism: Existing methods generally adopt a step-by-step loosely coupled solution approach. Taking the aforementioned method for predicting icing of rotating blades in aero-engines as an example, it also uses a sequential solution method to calculate the flow field, water droplet impact, and icing thermodynamics. It does not synchronously feed back the mass change and latent heat release caused by freezing phase change to the water film control equation for iteration. The lack of real-time two-way strong coupling between water film flow, heat transfer, and phase change leads to distortion of the coupling relationship of physical processes.

[0006] 3. Insufficient adaptation to the rotational conditions and three-dimensional surface features of wind turbines: Existing models are mostly designed for fixed-wing aircraft. For example, CN120724618B is mainly for rotating components of aero engines. Its control equations and coordinate system are not specifically designed for the large curvature three-dimensional surface of wind turbine blades. The influence of centrifugal force and Coriolis force caused by blade rotation on water film flow is not fully considered, and it cannot accurately reflect the unique water film distribution and freezing pattern of wind turbine blades.

[0007] 4. It is difficult to balance computational efficiency and accuracy: High-precision three-dimensional icing simulation requires a large amount of computation. Existing methods lack optimization in terms of time step control and mesh update strategies, making it difficult to meet the needs of engineering applications.

[0008] In summary, it is currently impossible to achieve a comprehensive solution that integrates precise simulation of three-dimensional curved surface water films for wind turbine blades, strong bidirectional coupling of flow-heat transfer-phase change, complete consideration of rotational effects, and a balance between accuracy and efficiency. Therefore, a high-precision water film flow-heat transfer-phase change coupled solution method that can overcome the above-mentioned shortcomings is urgently needed. Summary of the Invention

[0009] To address the problems of simplified water film models, loose coupling mechanisms, insufficient adaptation of rotational effects to three-dimensional curved surfaces, and difficulty in balancing accuracy and efficiency in existing technologies, this invention aims to provide a coupled solution method and system for water film flow-heat transfer-phase change during the icing process of wind turbine blades. By constructing a complete set of water film control equations in a three-dimensional curved surface body-fitted coordinate system, introducing rotational effects and realizing bidirectional strong coupling feedback between freezing phase change and control equations, and combining partitioned adaptive time step and dynamic mesh update, the invention improves computational efficiency while ensuring the accuracy of the physical mechanism, thereby achieving a more realistic, stable, and efficient numerical simulation of the icing process of wind turbine blades.

[0010] To achieve the above-mentioned technical objectives, the technical solution adopted by the present invention is as follows:

[0011] In a first aspect, the present invention discloses a coupled solution method for water film flow-heat transfer-phase change during the icing process of wind turbine blades, the method comprising the following steps:

[0012] S1: Generate a body-fitted computational mesh on the surface of the target wind turbine blade and establish a three-dimensional curved surface coordinate system, in which the mesh cells cover the pressure surface, suction surface and leading edge region of the blade.

[0013] S2: Obtain the flow field parameters and meteorological parameters at the current moment, and calculate the mass flow rate of water droplet impact on the blade surface based on the flow field parameters; where the flow field parameters include the pressure distribution on the blade surface, the airflow shear force and the convective heat transfer coefficient, and the meteorological parameters include the incoming flow temperature and the incoming flow liquid water content;

[0014] S3: Using the flow field parameters, meteorological parameters, and water droplet impact mass flow rate obtained in step S2 as input boundary conditions, solve the water film flow-heat transfer-phase change coupled control equations simultaneously on each grid cell. The equations include the water film mass conservation equation, the water film momentum conservation equation, and the water film energy conservation equation. The water film mass conservation equation is constructed with the water droplet impact mass as the source term and the evaporation or sublimation mass and the freezing phase change mass as the sink term. The water film momentum conservation equation is constructed based on the tangential forces on the three-dimensional curved surface, introducing airflow shear force, gravity component, centrifugal force component, Coriolis force component, and wall viscous drag. The water film energy conservation equation is constructed based on the energy balance relationship, introducing convective heat transfer, evaporative heat loss, latent heat of phase change, and external heat source.

[0015] S4: Based on the thermal balance calculation of the grid cell, the freezing ratio coefficient is determined, the freezing mass flow rate of the water film in each grid cell is determined, and the freezing mass flow rate is fed back to the water film mass conservation equation as a mass sink term. At the same time, the latent heat of phase change is fed back to the water film energy conservation equation as an energy source term to achieve a two-way strongly coupled iterative solution.

[0016] S5: Calculate the ice layer growth thickness of each grid cell within the current time step based on the freezing mass flow rate, and dynamically update the local geometry of the blade surface based on the ice layer growth thickness.

[0017] S6: Combining the partitioned adaptive time step strategy, based on the difference in icing rate in different regions of the blade, the first time step is used for the leading edge region, and the second time step is used for the pressure surface and suction surface. A smooth transition of time step is performed at the region boundary, wherein the first time step is smaller than the second time step. Steps S2 to S5 are repeated until the preset total icing time is reached, and the three-dimensional ice shape distribution and water film dynamic evolution process on the blade surface are output.

[0018] Furthermore, in step S3, the expression for the water film momentum conservation equation in the three-dimensional curved surface coordinate system is:

[0019] ;

[0020] Where t is time, The density of water, For water film thickness, The average velocity vector of the water film. Represents the vector dyadic product. For surface divergence operators, For airflow shear force, Let be the component of gravity tangentially to the curved surface. Let be the component of centrifugal acceleration in the tangential direction of the surface. For Coriolis force, The first term represents the dynamic viscosity of water, and the last term represents the wall friction resistance.

[0021] Furthermore, in step S3, the energy conservation equation for the water film is:

[0022] ;

[0023] Where t is time, The density of water, For water film thickness, The average velocity vector of the water film. For surface divergence operators, The specific heat capacity of water, The water film temperature, For convective heat transfer, heat flux density For evaporation or sublimation heat flux density, To freeze the latent heat flux density of the phase change, This represents the heat flux density of the external heat source.

[0024] Furthermore, in step S4, the freezing ratio coefficient is determined by solving the heat balance equation of the mesh element. :

[0025] ;

[0026] in, The specific heat capacity of water, The water film temperature, For convective heat transfer, heat flux density For evaporation or sublimation heat flux density, The heat flux density of the external heat source. The mass flow rate of water droplet impact per unit area. This is the freezing temperature. For the latent heat of freezing of water; when At that time, take , indicates no freeze; when At that time, take , indicates that everything is frozen.

[0027] Further, in step S5, the ice layer thickness growth is calculated using the following formula. :

[0028] ;

[0029] in, The density of ice, For the current time step, This is the freezing ratio coefficient. The mass flow rate is the impact mass of water droplets per unit area.

[0030] Step S6 further includes:

[0031] S61: Divide the blade surface into multiple characteristic regions, including at least the leading edge stagnation region, the pressure surface laminar flow region, and the suction surface separation region;

[0032] S62: Based on the differences in icing rate of each feature region, set differentiated initial time steps for different feature regions, with the leading edge stagnation region using the minimum initial time step;

[0033] S63: After the solution is completed at each time step, monitor the equation residuals and physical quantity gradients of each feature region. When the gradient of any feature region exceeds the preset gradient threshold, automatically reduce the time step size of the next iteration for that feature region.

[0034] S64: When the ratio of time steps of adjacent feature regions exceeds a preset ratio threshold, the time step is linearly smoothed at the region boundary to ensure the stability of the numerical solution.

[0035] Furthermore, in step S6, before each iteration, the three-dimensional curved surface body mesh is reconstructed based on the updated ice layer thickness and the local geometry of the blade surface, and the external flow field is resolved to update the flow field parameters and boundary conditions on the blade surface, thereby achieving fully coupled iteration of the evolution of the flow field, water film and ice shape.

[0036] Secondly, this invention discloses a coupled solution system for water film flow-heat transfer-phase change during the icing process of wind turbine blades, the system comprising:

[0037] The mesh generation module is used to generate a body-fitted computational mesh on the surface of the target wind turbine blade and establish a three-dimensional curved surface coordinate system, in which the mesh cells cover the blade pressure surface, suction surface and leading edge region;

[0038] The data acquisition interface is used to acquire the flow field parameters and meteorological parameters at the current moment, and calculate the mass flow rate of water droplet impact on the blade surface based on the flow field parameters; among which, the flow field parameters include the pressure distribution on the blade surface, the airflow shear force and the convective heat transfer coefficient, and the meteorological parameters include the incoming flow temperature and the incoming flow liquid water content;

[0039] The core module for solving the water film problem uses the acquired flow field parameters, meteorological parameters, and water droplet impact mass flow rate as input boundary conditions to solve the coupled control equations of water film flow-heat transfer-phase change on each grid cell. The equations include the water film mass conservation equation, the water film momentum conservation equation, and the water film energy conservation equation. The water film mass conservation equation is constructed with the water droplet impact mass as the source term and the evaporation or sublimation mass and the freezing phase change mass as the sink term. The water film momentum conservation equation is constructed based on the tangential forces on the three-dimensional curved surface, introducing airflow shear force, gravity component, centrifugal force component, Coriolis force component, and wall viscous drag. The water film energy conservation equation is constructed based on the energy balance relationship, introducing convective heat transfer, evaporative heat loss, latent heat of phase change, and external heat sources.

[0040] The freezing calculation module is used to calculate the freezing ratio coefficient based on the thermal balance of grid cells, determine the freezing mass flow rate of the water film in each grid cell, calculate the ice layer growth thickness of each grid cell within the current time step based on the freezing mass flow rate, and feed the freezing mass flow rate as a mass sink term back to the water film mass conservation equation, while feeding the latent heat of phase change as an energy source term back to the water film energy conservation equation, so as to achieve a two-way strongly coupled iterative solution.

[0041] The mesh update module is used to dynamically update the local geometry of the blade surface based on the ice thickness.

[0042] The time step control module is used to combine the partitioned adaptive time step strategy, and according to the difference in icing rate in different regions of the blade, the first time step is used for the leading edge region, the second time step is used for the pressure surface and the suction surface, and the time step is smoothly transitioned at the region boundary, wherein the first time step is smaller than the second time step.

[0043] The data output module is used to output the three-dimensional ice shape distribution on the blade surface and the dynamic evolution process of the water film, including the water film thickness field, temperature field, and flow velocity field.

[0044] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0045] First, the water film flow-heat transfer-phase change coupled solution method and system of the present invention in the wind turbine blade icing process fully describes the two-way coupling mechanism among water film flow, heat transfer and phase change, and simultaneously feeds back the freezing ratio coefficient to the mass equation and energy equation, realizing the two-way strong coupling of water film flow-heat transfer-phase change, with a more complete physical mechanism, and can accurately predict complex icing phenomena such as overflow water and ice corner formation. It is the first time that a three-dimensional curved surface water film flow-heat transfer-phase change fully coupled solution has been realized for wind turbine blades.

[0046] Secondly, the water film flow-heat transfer-phase change coupled solution method and system of the present invention in the wind turbine blade icing process fully considers centrifugal force and Coriolis force in the momentum equation, and can accurately simulate the influence of wind turbine rotation on water film distribution. This is a key innovation that distinguishes it from aviation icing models.

[0047] Third, the water film flow-heat transfer-phase change coupled solution method and system of the present invention in the process of wind turbine blade icing improves the overall calculation efficiency by more than 50% while ensuring the calculation accuracy of the critical leading edge region through a differentiated time step strategy. The partitioned adaptive time step strategy adopted in the present invention can significantly improve efficiency.

[0048] Fourth, the water film flow-heat transfer-phase change coupled solution method and system of the present invention for the wind turbine blade icing process adopts dynamic mesh update in real time as the ice layer grows to ensure long-term simulation accuracy, avoiding the error accumulation caused by excessive geometric deformation in traditional methods.

[0049] Fifth, the water film flow-heat transfer-phase change coupled solution method and system of the present invention for the wind turbine blade icing process can be integrated into the existing CFD framework, providing accurate input for anti-icing system design, ice shape database construction, and icing detector layout optimization, and has strong engineering applicability. Attached Figure Description

[0050] Figure 1 This is an overall flowchart of the water film flow-heat transfer-phase change coupled solution method for the wind turbine blade icing process of the present invention.

[0051] Figure 2 This is a dynamic mesh update diagram of the airfoil section after the blades have iced.

[0052] Figure 3 This is a schematic diagram illustrating the solution logic for the water film governing equations.

[0053] Figure 4 This is an example diagram of three-dimensional icing of a blade simulated using the method of the present invention. Detailed Implementation

[0054] The embodiments of the present invention will be described in further detail below with reference to the accompanying drawings.

[0055] See Figure 1 This invention discloses a coupled solution method for water film flow-heat transfer-phase change during the icing process of wind turbine blades, the method comprising the following steps:

[0056] S1: Generate a body-fitted computational mesh on the surface of the target wind turbine blade and establish a three-dimensional curved surface coordinate system, in which the mesh cells cover the pressure surface, suction surface and leading edge region of the blade.

[0057] S2: Obtain the flow field parameters and meteorological parameters at the current moment, and calculate the mass flow rate of water droplet impact on the blade surface based on the flow field parameters; where the flow field parameters include the pressure distribution on the blade surface, the airflow shear force and the convective heat transfer coefficient, and the meteorological parameters include the incoming flow temperature and the incoming flow liquid water content;

[0058] S3: Using the flow field parameters, meteorological parameters, and water droplet impact mass flow rate obtained in step S2 as input boundary conditions, solve the water film flow-heat transfer-phase change coupled control equations simultaneously on each grid cell. The equations include the water film mass conservation equation, the water film momentum conservation equation, and the water film energy conservation equation. The water film mass conservation equation is constructed with the water droplet impact mass as the source term and the evaporation or sublimation mass and the freezing phase change mass as the sink term. The water film momentum conservation equation is constructed based on the tangential forces on the three-dimensional curved surface, introducing airflow shear force, gravity component, centrifugal force component, Coriolis force component, and wall viscous drag. The water film energy conservation equation is constructed based on the energy balance relationship, introducing convective heat transfer, evaporative heat loss, latent heat of phase change, and external heat source.

[0059] S4: Based on the thermal balance calculation of the grid cell, the freezing ratio coefficient is determined, the freezing mass flow rate of the water film in each grid cell is determined, and the freezing mass flow rate is fed back to the water film mass conservation equation as a mass sink term. At the same time, the latent heat of phase change is fed back to the water film energy conservation equation as an energy source term to achieve a two-way strongly coupled iterative solution.

[0060] S5: Calculate the ice layer growth thickness of each grid cell within the current time step based on the freezing mass flow rate, and dynamically update the local geometry of the blade surface based on the ice layer growth thickness.

[0061] S6: Combining the partitioned adaptive time step strategy, according to the difference in icing rate in different regions of the blade, set a time step that matches the region, and perform a smooth transition of time step at the region boundary; repeat steps S2 to S5 until the preset total icing time is reached, and output the three-dimensional ice shape distribution and water film dynamic evolution process on the blade surface.

[0062] The proposed method for coupled solution of water film flow, heat transfer and phase change during the icing process of wind turbine blades, with "establishment of water film control equations - fully coupled solution - adaptive time advancement - geometric dynamic update" as the main line, has for the first time achieved accurate simulation of three-dimensional curved surface water film for wind turbine blades.

[0063] Step 1: Body-fit computational mesh generation and coordinate system establishment

[0064] A structured or unstructured body-fitted computational mesh is generated on the surface of the target blade, and a three-dimensional surface coordinate system is established. ,in , These are the tangential coordinates of the surface (corresponding to the chord and span directions of the blade, respectively). The coordinates are normal coordinates. Water film thickness. along Directional measurement, average velocity of water film for , The component vector of direction.

[0065] Step 2: Establishing the Coupled Control Equations

[0066] For each grid cell, establish the following set of governing equations:

[0067] ① Water film mass conservation equation:

[0068] ;

[0069] in, The mass flow rate of water droplet impact per unit area is used as the water droplet impact mass source term (provided by the water droplet trajectory calculation module). This is the mass summation of evaporation / sublimation (determined by heat transfer calculations). This is a frozen quality pool (to be solved).

[0070] ② Equation for conservation of momentum of water film:

[0071] The equation for the conservation of water film momentum in a three-dimensional surface coordinate system is expressed as follows:

[0072] ;

[0073] in, The density of water, For water film thickness, The average velocity vector of the water film. For surface divergence operators, For airflow shear force, Let be the component of gravity tangentially to the curved surface. Let be the component of centrifugal acceleration in the tangential direction of the surface. For Coriolis force, The first term represents the dynamic viscosity of water, and the last term represents the wall friction resistance.

[0074] The momentum conservation equation for water film considers airflow shear force (the main driving force for water film flow), the tangential component of gravity, the tangential component of centrifugal force, Coriolis force, and wall friction resistance. Centrifugal acceleration... ,in The blade rotational angular velocity vector. This is the position vector from the center of the grid cell to the rotation axis.

[0075] ③ Water film energy conservation equation:

[0076] The energy conservation equation for water film is:

[0077] ;

[0078] in, The specific heat capacity of water, The water film temperature, For convective heat transfer, heat flux density Evaporation / sublimation heat flux density, The heat flux density for the latent heat of the freezing phase change (a positive value indicates the release of latent heat). Heat flux density from an external heat source (such as electric heating)

[0079] The water film energy conservation equation of this invention considers convective heat transfer (coupled with the flow field), evaporative heat loss, latent heat of freezing phase change, and external heat sources. Among these, the latent heat of freezing phase change heat flux density... This reflects the feedback of the freezing process on the energy balance.

[0080] Step 3: Solving the freezing ratio coefficient and bidirectional coupling

[0081] Solving for the freezing ratio coefficient based on the thermal equilibrium state of the mesh cells. :

[0082] ;

[0083] in, The mass flow rate of water droplet impact per unit area. This is the freezing temperature (usually taken as 0℃). For the latent heat of freezing of water; when At that time, take , indicates no freeze; when At that time, take This indicates a complete freeze. Freeze percentage coefficient. The physical meaning is: what percentage of the water droplets impacting this unit freezes immediately, while the remaining portion continues to flow as liquid water. This invention's freezing ratio coefficient... The calculations rely on both the energy equation (via the heat flow term) and the mass equation (via the heat flow term). ), and the amount of frozen It is also fed back as a source term to the mass and energy equations, thus achieving bidirectional coupling.

[0084] Step 4: Calculation of ice layer growth thickness and geometric update

[0085] The ice layer thickness is calculated based on the freezing mass. The calculation formula is:

[0086] ;

[0087] in, The density of ice, This is the current time step. Then, based on the ice layer thickness increase... Update the geometry at the grid cell. When the ice layer thickness accumulates to a certain level (e.g., exceeding 20% ​​of the grid cell height), trigger the grid update module to reconstruct the local mesh and resolve the flow field.

[0088] Step 5: Partition Adaptive Time Step Advancement

[0089] To address the timescale differences in the icing process across different regions of wind turbine blades, a zoned adaptive time step strategy is adopted. Based on the varying icing rates across different regions of the blade, a matching time step is set for each region. Specifically:

[0090] S61: Divide the blade surface into multiple characteristic regions, including at least the leading edge stagnation region, the pressure surface laminar flow region, and the suction surface separation region;

[0091] S62: Based on the differences in icing rate of each feature region, set differentiated initial time steps for different feature regions, with the leading edge stagnation region using the minimum initial time step;

[0092] S63: After the solution is completed at each time step, monitor the equation residuals and physical quantity gradients of each feature region. When the gradient of any feature region exceeds the preset gradient threshold, automatically reduce the time step size of the next iteration for that feature region.

[0093] S64: When the ratio of time steps of adjacent feature regions exceeds a preset ratio threshold, the time step is linearly smoothed at the region boundary to ensure the stability of the numerical solution.

[0094] Combined with appendix Figures 1-4 The solution method and the matching coupled solution system of the present invention will be described in detail. Figure 1 This is a flowchart illustrating the overall process of the coupled solution in this invention. Figure 2 This is a schematic diagram of the dynamic mesh update of the airfoil section after icing. Figure 3 This is a schematic diagram illustrating the logic for solving the water film governing equations. Figure 4 This is a simulation image of the three-dimensional icing morphology of the blade obtained by this method.

[0095] All embodiments below use a fixed physical property parameter: water density. =1000 ice density =917 The specific heat capacity of liquid water at constant pressure is 4200. Latent heat of water upon freezing = 335 Freezing temperature of ice-water phase transition =0℃, and the physical property parameters will not be repeated thereafter.

[0096] Example 1: Three-dimensional icing simulation under typical working conditions

[0097] This embodiment takes a 2MW wind turbine blade as the object and simulates the formation, migration, freezing and ice layer growth process of water film on the blade surface under typical low temperature and humid inflow conditions to illustrate the specific implementation process of the method of the present invention.

[0098] Step 1: Establishing the blade model and computational grid

[0099] A three-dimensional geometric model of the wind turbine blade to be analyzed was selected as the computational object. This blade exhibits variable chord length and variable twist angle along its span, with a chord length of approximately 2.5 m near the blade root and approximately 1.2 m near the blade tip. The blade twist angle gradually changes from the blade root to the blade tip. Based on the blade geometric model, a body-fitted computational mesh was generated on the blade surface, covering the leading edge region, pressure surface region, suction surface region, and the high-curvature region near the blade tip. To improve the computational resolution of the leading edge icing region, localized mesh refinement was applied near the leading edge; a relatively coarser mesh scale was used for regions far from the leading edge with less water droplet impact. The total number of meshes in the computational domain is approximately 2 million, of which approximately 50,000 are on the blade surface.

[0100] A three-dimensional surface coordinate system is established on the blade surface, with the chordal and spanwise directions used as the tangential coordinates and the direction perpendicular to the blade surface used as the normal coordinates. The water film thickness is measured along the normal direction, and the water film velocity is decomposed into chordal and spanwise components. This surface coordinate system allows for the discrete solution of subsequent equations for water film mass conservation, momentum conservation, and energy conservation on the complex surface of the blade.

[0101] Step 2: Obtaining flow field parameters and water droplet impact parameters

[0102] Environmental conditions were set as follows: incoming wind speed 10 m / s, ambient temperature -5℃, liquid water content in the air 0.3 g / m³, median water droplet diameter 20 μm, and rated blade rotation speed 15 rpm. The SST k-ω turbulence model was used to solve the RANS external flow field, extracting the global surface pressure, airflow shear force, and convective heat transfer coefficient of the blades. The local water droplet collection efficiency β was calculated based on the water droplet trajectory algorithm. Solve for the mass flow rate of water droplet impact in each grid cell. The mass flow rate of water droplet impact per unit area. The free-flow velocity is at infinity. Simulation results show that the water droplet impact flow rate in the stagnation zone at the leading edge of the blade is much higher than that in the far regions of the pressure and suction surfaces.

[0103] Step 3: Solving the water film-heat transfer-phase change coupling

[0104] Using flow field parameters, meteorological parameters, and water droplet impact flow rate as boundary conditions, the aforementioned water film mass, momentum, and energy conservation control equations are simultaneously solved discretly; the freezing ratio coefficient of each grid cell is calculated based on the heat balance formula. The mass flow rate generated during freezing is fed back to the water film mass conservation equation as a sink term, and the latent heat of phase change is fed back to the water film energy conservation equation as an energy source term, thereby achieving strong bidirectional coupling of flow, heat transfer, and phase change.

[0105] Taking a single mesh cell at the front edge as an example, the mesh freezing ratio at 120s during simulation. =0.65, 65% of the impacting water droplets freeze in situ, and the remaining liquid water is transported to the upstream and downstream grids under the combined action of airflow shear force, rotational centrifugal force, and Coriolis force, realizing the simulation of overflow water migration and refreezing.

[0106] Step 4: Ice layer growth and local geometric update

[0107] The ice thickness increase per grid step was calculated using the ice thickness calculation formula. The local geometry is updated along the blade normal. In this embodiment, the reference height of the grid normal is set to 0.2 mm. When the ice layer thickening reaches 20% of the grid height (critical value 0.04 mm), local body-fitting grid reconstruction is automatically triggered, and the boundary parameters of the surrounding flow field are corrected simultaneously. The dynamic grid change shape is as follows: Figure 2 As shown.

[0108] Step 5: Results Output and Analysis

[0109] After a total simulation time of 1800 seconds, the three-dimensional ice shape, water film thickness field, temperature field, and velocity field data of the entire blade area were output. The ice layer was mainly concentrated at the leading edge of the blade. Under the influence of centrifugal force, liquid water overflowed towards the blade tip along the spanwise direction and locally refrozen to form ice corners. The final three-dimensional ice shape simulation results are shown below. Figure 4 The blue area represents the ice layer on the leaf surface. Refer to the appendix for the iterative logic of the entire equation system. Figure 3 .

[0110] Example 2: Analysis of the Influence of Rotation Effect on Water Film Migration and Ice Shape Distribution

[0111] Using the same blade geometry and mesh scheme as in Example 1, the operating conditions were set as follows: incoming air velocity 8 m / s, incoming air temperature -8℃, liquid water content 0.4 g / m³, blade rotation speed 18 rpm, and total icing simulation time of 1800 s. Two control schemes, A and B, were set up.

[0112] Group A adopts the water film momentum conservation equation described in this invention, and considers airflow shear force, gravity tangential component, centrifugal force component, Coriolis force component and wall resistance in the force term of the water film;

[0113] Group B, while keeping other conditions constant, does not consider the centrifugal force component and Coriolis force component caused by blade rotation, but only the airflow shear force, the tangential component of gravity, and the wall drag.

[0114] Simulation data comparison: Group A takes into account rotational inertia force, and centrifugal force drives the water film to flow towards the blade tip along the span. After the simulation, the average ice layer thickness at the blade tip is 12.3 mm and the ice layer thickness at the blade root is 6.1 mm, with a significant ice thickness gradient along the span. Group B has no rotational inertia force, and the water film migration along the span is extremely weak. The ice layer thickness range of the entire blade is only 7.2 to 7.8 mm, and the ice thickness is evenly distributed along the span.

[0115] The comparative results demonstrate that introducing centrifugal force and Coriolis force into the simulation of icing on wind turbine rotating blades can accurately reproduce the water film overflow pattern under rotating conditions, thereby improving the realism of the ice shape simulation.

[0116] Example 3: Calculation process of the partitioned adaptive time step strategy

[0117] Using the blade model and meteorological conditions from Example 1, the blade was divided into three major characteristic zones: the leading edge stagnation zone, the pressure surface region, and the suction surface region. Comparative tests were conducted on three types of time-step control strategies, with pre-set quantitative control thresholds: water film thickness gradient threshold. The ratio of time steps between adjacent partitions is capped at 3 times. If the threshold is exceeded, the step size will be automatically adjusted or the boundary will be smoothly interpolated.

[0118] To verify the effectiveness of the strategy, three time-stepping methods were compared: Strategy 1 used a globally fixed small time step Δt=0.01s, which had the best calculation accuracy, but the overall simulation time was 12.6h; Strategy 2 used a globally fixed large time step Δt=0.05s, which took 2.1h, but the simulation error of the leading edge ice shape was 19.7%, and the accuracy did not meet the engineering requirements; Strategy 3 adopted the partitioned adaptive time stepping strategy of this invention, with an initial stepping of 0.01s for the leading edge and 0.05s for the pressure and suction surfaces; after each step calculation, the gradient of physical quantities in the partition and the iteration residual were monitored, and the time stepping of the partition that exceeded the limit was automatically reduced; when the stepping ratio of adjacent partitions exceeded 3 times the threshold, a linear smooth transition was achieved at the boundary.

[0119] The final calculation time of the proposed solution was 5.8 hours, which is 53.9% more efficient than the global small step size solution. The overall ice shape simulation error was less than 3%. The calculation cycle was significantly shortened while ensuring the calculation accuracy of key icing areas, achieving a balance between accuracy and efficiency.

[0120] Example 4: Implementation of a Coupled Solution System for Wind Turbine Blade Icing

[0121] The solution system of this invention is based on secondary development of general-purpose CFD software and deployed on a 16-core industrial computing workstation, realizing fully automated calculation without human intervention:

[0122] 1. Mesh generation module: Interfaces with ICEM and OpenFOAM meshes to automatically generate body-fitted meshes for blade surfaces and automatically refine the mesh in the leading edge region;

[0123] 2. Data acquisition interface: Connects to an external flow field solver and meteorological parameter input port, automatically imports wind speed, temperature, and liquid water content, and calculates the mass flow rate of water droplet impact in the entire area in batches;

[0124] 3. Water film solution core module + freezing calculation module: It has three built-in subroutines for discrete solution of conservation equations and code for solving freezing ratio coefficients, and automatically completes bidirectional feedback iteration of freezing quality and latent heat of phase change;

[0125] 4. Mesh Update Module: Embedded dynamic mesh deformation and local reconstruction scripts automatically update blade geometry and surface mesh after reaching the ice thickness threshold;

[0126] 5. Time step control module: integrates a partitioned adaptive step size algorithm, which automatically adjusts the calculation step size of each region and performs boundary smooth interpolation according to the partition threshold;

[0127] 6. Data Output Module: Outputs 3D ice shape, water film thickness / temperature / flow velocity field data in Tecplot format, which can be directly used for post-processing plotting.

[0128] In actual operation of Example 1, the system automatically completes the process from mesh generation to result output without manual intervention, demonstrating strong engineering adaptability.

[0129] Industrial Application Example: Optimized Design of Electrothermal Anti-icing and De-icing System for Wind Turbine Blades

[0130] A domestic wind power equipment manufacturer used the solution system of this invention to optimize the selection of power for electrothermal anti-icing of new blades. With the meteorological environment kept constant, four levels of heating power density were set for simulation comparison: 0 (no heating), 2kW / m², 4kW / m², and 6kW / m², with a single simulation duration of 30 minutes.

[0131] • No heating: The average ice thickness at the leading edge of the blade is 15mm, resulting in sharp ice angles;

[0132] • 2kW / m²: The average ice thickness is 8mm, and local ice corner defects still exist;

[0133] • 4kW / m²: The average ice layer thickness is 3mm, with no ice corners formed and the surface water film flows smoothly;

[0134] ·6kW / m²: No ice formation on the blade surface, but high heating energy consumption and excessive operating costs.

[0135] Based on the simulation results, the manufacturer selected 4kW / m² as the rated heating power for the mass-produced model and conducted a field low-temperature ice-attaching test on the prototype. The prototype was continuously run for 30 minutes in an actual environment of -5℃ and liquid water content of 0.3g / m³. The actual ice layer thickness on the blades was 2.7-3.2mm, and the error between the simulation results and the actual results was less than 10%, thus verifying the reliability of the simulation data through actual testing.

[0136] Although preferred embodiments of this application have been described, those skilled in the art, upon learning the basic inventive concept, can make other changes and modifications to these embodiments. Therefore, the appended claims are intended to be interpreted as including the preferred embodiments as well as all changes and modifications falling within the scope of this application.

[0137] Obviously, those skilled in the art can make various modifications and variations to this application without departing from the spirit and scope of this application. Therefore, if such modifications and variations fall within the scope of the claims of this application and their equivalents, this application also intends to include such modifications and variations.

Claims

1. A coupled solution method for water film flow-heat transfer-phase change during the icing process of wind turbine blades, characterized in that, The method includes the following steps: S1: Generate a body-fitted computational mesh on the surface of the target wind turbine blade and establish a three-dimensional curved surface coordinate system, in which the mesh cells cover the pressure surface, suction surface and leading edge region of the blade. S2: Obtain the flow field parameters and meteorological parameters at the current moment, and calculate the mass flow rate of water droplet impact on the blade surface based on the flow field parameters; where the flow field parameters include the pressure distribution on the blade surface, the airflow shear force and the convective heat transfer coefficient, and the meteorological parameters include the incoming flow temperature and the incoming flow liquid water content; S3: Using the flow field parameters, meteorological parameters, and water droplet impact mass flow rate obtained in step S2 as input boundary conditions, solve the water film flow-heat transfer-phase change coupled control equations simultaneously on each grid cell. The equations include the water film mass conservation equation, the water film momentum conservation equation, and the water film energy conservation equation. The water film mass conservation equation is constructed with the water droplet impact mass as the source term and the evaporation or sublimation mass and the freezing phase change mass as the sink term. The water film momentum conservation equation is constructed based on the tangential forces on the three-dimensional curved surface, introducing airflow shear force, gravity component, centrifugal force component, Coriolis force component, and wall viscous drag. The water film energy conservation equation is constructed based on the energy balance relationship, introducing convective heat transfer, evaporative heat loss, latent heat of phase change, and external heat source. S4: Based on the thermal balance calculation of the grid cell, the freezing ratio coefficient is determined, the freezing mass flow rate of the water film in each grid cell is determined, and the freezing mass flow rate is fed back to the water film mass conservation equation as a mass sink term. At the same time, the latent heat of phase change is fed back to the water film energy conservation equation as an energy source term to achieve a two-way strongly coupled iterative solution. S5: Calculate the ice layer growth thickness of each grid cell within the current time step based on the freezing mass flow rate, and dynamically update the local geometry of the blade surface based on the ice layer growth thickness. S6: Combining the partitioned adaptive time step strategy, the blade surface is divided into multiple feature regions, including at least the leading edge stagnation region, the pressure surface region, and the suction surface region; according to the icing rate, physical quantity gradient, or equation residual of each region, the corresponding local time step is set for different regions, and the time step is smoothly transitioned between adjacent regions; repeat steps S2 to S5 until the preset total icing time is reached, and output the three-dimensional ice shape distribution and water film dynamic evolution process of the blade surface.

2. The coupled solution method for water film flow-heat transfer-phase change during the icing process of wind turbine blades according to claim 1, characterized in that, In step S3, the expression for the water film momentum conservation equation in the three-dimensional curved surface coordinate system is: ； Where t is time, The density of water, For water film thickness, The average velocity vector of the water film. Represents the vector dyadic product. For surface divergence operators, For airflow shear force, Let be the component of gravity tangentially to the curved surface. Let be the component of centrifugal acceleration in the tangential direction of the surface. For Coriolis force, The first term represents the dynamic viscosity of water, and the last term represents the wall friction resistance.

3. The coupled solution method for water film flow-heat transfer-phase change during the icing process of wind turbine blades according to claim 1, characterized in that, In step S3, the energy conservation equation for the water film is: ； Where t is time, The density of water, For water film thickness, The average velocity vector of the water film. For surface divergence operators, The specific heat capacity of water, The water film temperature, For convective heat transfer, heat flux density For evaporation or sublimation heat flux density, To freeze the latent heat flux density of the phase change, This represents the heat flux density of the external heat source.

4. The coupled solution method for water film flow-heat transfer-phase change during the icing process of wind turbine blades according to claim 1, characterized in that, In step S4, the freezing ratio coefficient is determined by solving the heat balance equation of the mesh element. : ； in, The specific heat capacity of water, The water film temperature, For convective heat transfer, heat flux density For evaporation or sublimation heat flux density, The heat flux density of the external heat source. The mass flow rate of water droplet impact per unit area. This is the freezing temperature. For the latent heat of freezing of water; when At that time, take , indicates no freeze; when At that time, take , indicates that everything is frozen.

5. The coupled solution method for water film flow-heat transfer-phase change during the icing process of wind turbine blades according to claim 1, characterized in that, In step S5, the ice layer thickness growth is calculated using the following formula. : ； in, The density of ice, For the current time step, This is the freezing ratio coefficient. The mass flow rate is the water droplet impact mass per unit area.

6. The coupled solution method for water film flow-heat transfer-phase change during the icing process of wind turbine blades according to claim 1, characterized in that, Step S6 further includes: S61: Divide the blade surface into multiple characteristic regions, including at least the leading edge stagnation region, the pressure surface laminar flow region, and the suction surface separation region; S62: Based on the differences in icing rate of each feature region, set differentiated initial time steps for different feature regions, with the leading edge stagnation region using the minimum initial time step; S63: After the solution is completed at each time step, monitor the equation residuals and physical quantity gradients of each feature region. When the gradient of any feature region exceeds the preset gradient threshold, automatically reduce the time step size of the next iteration for that feature region. S64: When the ratio of time steps of adjacent feature regions exceeds a preset ratio threshold, the time step is linearly smoothed at the region boundary to ensure the stability of the numerical solution.

7. The coupled solution method for water film flow-heat transfer-phase change during the icing process of wind turbine blades according to claim 1, characterized in that, In step S6, when the ice layer thickness increment, local curvature change, or mesh quality index meets the preset update conditions, the three-dimensional curved surface body mesh is reconstructed based on the updated ice layer thickness and the local geometry of the blade surface, and the external flow field is resolved or locally modified to update the flow field parameters and boundary conditions on the blade surface, thereby realizing the fully coupled iteration of flow field, water film, and ice shape evolution.

8. A coupled solution system for water film flow-heat transfer-phase change during the icing process of wind turbine blades based on the method of any one of claims 1-7, characterized in that, The system includes: The mesh generation module is used to generate a body-fitted computational mesh on the surface of the target wind turbine blade and establish a three-dimensional curved surface coordinate system, in which the mesh cells cover the blade pressure surface, suction surface and leading edge region. The data acquisition interface is used to acquire the flow field parameters and meteorological parameters at the current moment, and calculate the mass flow rate of water droplet impact on the blade surface based on the flow field parameters; among which, the flow field parameters include the pressure distribution on the blade surface, the airflow shear force and the convective heat transfer coefficient, and the meteorological parameters include the incoming flow temperature and the incoming flow liquid water content; The core module for solving the water film problem uses the acquired flow field parameters, meteorological parameters, and water droplet impact mass flow rate as input boundary conditions. It simultaneously solves the water film flow-heat transfer-phase change coupled control equations on each grid cell. The equations include the water film mass conservation equation, the water film momentum conservation equation, and the water film energy conservation equation. The water film mass conservation equation is constructed with the water droplet impact mass as the source term and the evaporation or sublimation mass and the freezing phase change mass as the sink term. The water film momentum conservation equation is constructed based on the tangential forces on the three-dimensional curved surface, introducing airflow shear force, gravity component, centrifugal force component, Coriolis force component, and wall viscous drag. The water film energy conservation equation is constructed based on the energy balance relationship, introducing convective heat transfer, evaporative heat loss, latent heat of phase change, and external heat sources. The freezing calculation module is used to calculate the freezing ratio coefficient based on the thermal balance of grid cells, determine the freezing mass flow rate of the water film in each grid cell, calculate the ice layer growth thickness of each grid cell within the current time step based on the freezing mass flow rate, and feed the freezing mass flow rate as a mass sink term back to the water film mass conservation equation, while feeding the latent heat of phase change as an energy source term back to the water film energy conservation equation, so as to achieve a two-way strongly coupled iterative solution. The mesh update module is used to dynamically update the local geometry of the blade surface based on the ice thickness. The time step control module is used to combine the partitioned adaptive time step strategy, set a matching time step for different regions according to the differences in icing rate in different regions of the blade, and perform a smooth transition of time step at the region boundary. The data output module is used to output the three-dimensional ice shape distribution on the blade surface and the dynamic evolution process of the water film, including the water film thickness field, temperature field, and flow velocity field.

Images