A method for predicting flutter of a turbine blade

Abstract

This invention discloses a method for predicting turbine blade flutter. The key point is the coupling of a real physical model of unequilibrium wet steam with the influence coefficient method for predicting flutter in the last-stage long blades of a turbine. This involves solid modeling of the last-stage moving blade, modal analysis, and harmonic response analysis, and constructing a nine-channel fluid domain for influence coefficient method calculations. The extracted main mode shapes are then interpolated onto the CFD nodes of the nine-channel fluid domain using a three-dimensional linear interpolation method. Furthermore, by modifying the wet steam model to accurately characterize the two-phase flow characteristics of the last-stage wet steam, unsteady flow field calculations using the nine-channel influence coefficient method are performed to obtain the aerodynamic damping coefficients for different nodal diameters. This enables the prediction of the flutter stability of the last-stage moving blades and the identification of critical nodal diameters. This significantly reduces the time cost of recalculating the tip phase angle for each blade throughout the entire revolution and improves the physical reliability and reproducibility of the prediction.

Description

Technical Field

[0001] This invention belongs to the field of flutter prediction technology, and in particular relates to a method for predicting flutter in steam turbine blades. Background Technology

[0002] In high-speed rotating machinery, the risk of flutter is particularly concentrated in the last stage blades of the low-pressure cylinder. The blades in this part are large in size, have high load, and operate in an extremely harsh environment: during the process of violent expansion, the steam crosses the saturation line and forms a typical non-equilibrium wet steam two-phase flow. This process is accompanied by irreversible phase change, latent heat release, condensation shock wave and droplet generation, which triggers a complex heat-fluid-structure coupling effect, which constitutes a continuous and strong unsteady disturbance to aerodynamic stability.

[0003] Current methods for predicting blade flutter mainly rely on indirect decoupling analysis methods such as Traveling Wave Mode (TWM). Existing technologies closely related to this invention cover two main areas: wet steam effect research and aeroelastic analysis methods. Regarding the wet steam effect, Fuhrer et al. (Fuhrer F, Seume JR The influence of non-equilibrium wet steam effects on the aeroelastic stability of steam turbineblades. ASME J Eng Gas Turbines Power, 2016, 138(6): 062601.) confirmed through comparative studies that the transient excitation caused by condensation shock waves under the non-equilibrium wet steam model has a significant impact on the local damping characteristics. Yao et al. (Yao Q, et al. Shock waves characteristics and losses estimation in wet steamflows. Int J Heat Mass Transf, 2025.) and Hu et al. (Hu Z, et al. Investigation of heterogeneous condensation and wave propagation. Energy Convers Manag, 2024.) further pointed out that the propagation behavior of condensation shock waves and the droplet size and distribution characteristics can significantly modulate the local pressure field, which can easily induce unsteady disturbances in the blade tip region.In terms of aeroelastic analysis methods, TWM assesses stability by calculating unsteady aerodynamic work, while ICM (Influence Coefficient Method) constructs the aerodynamic response under arbitrary traveling wave modes through linear superposition. Liu et al. (Liu Y., Chen C., Wang J. Verification of influence coefficient method with traveling wave mode under low interblade phaseangles. Journal of Thermal Science, 2020, 29(5): 933–941.) systematically compared the differences between ICM and TWM in aerodynamic damping calculation based on experimental data and numerical simulations. The results showed that the two methods have good consistency under most operating conditions. Phan and He (Phan TQ, He L. Comparison of traveling wavemode and influence coefficient methods in flutter prediction. ASME Turbo Expo2015, GT2015-42987.) further pointed out that the two methods differ in boundary condition handling, and ICM shows higher efficiency when multi-phase analysis is required.

[0004] The research objective of existing technologies is to accurately model and predict the flutter behavior of the last-stage blades under wet steam and variable operating conditions. The main technical measures are to use unsteady CFD (Computational Fluid Dynamics) simulation combined with a developed non-equilibrium wet steam model to capture complex flow fields, and to use the ICM indirect coupling method to calculate the aerodynamic damping under full pitch diameter. This coupling method can accurately reflect the true aeroelastic characteristics of the last-stage turbine under wet steam two-phase flow conditions.

[0005] However, existing technologies still have the following shortcomings:

[0006] 1) Disconnect between physical model and actual operating conditions: During the intense expansion process, the last-stage blades of a steam turbine are in a non-equilibrium wet steam two-phase flow environment. This process is accompanied by complex irreversible phase transitions, latent heat release, and condensation shock waves, resulting in strong thermal-fluid-structure coupled unsteady disturbances. However, traditional flutter prediction models still use ideal gas or other single-phase flows as the working fluid. This physical model, detached from actual operating conditions, completely ignores the precise influence of non-equilibrium condensation phase transitions on aerodynamic damping, leading to inherent biases in flutter prediction results and consequently failing to accurately assess the true aeroelastic stability of the last-stage blades under wet steam conditions.

[0007] 2) Limitations in prediction range and efficiency. Traditional TWM (Traveling Wave Mode) algorithms require different minimum flow paths for modeling different nodal diameter modes, some even requiring full-circle modeling, and each calculation can only obtain aerodynamic damping results for one nodal diameter. Due to limitations in computational cost and efficiency, many existing flutter prediction studies often only examine aerodynamic damping values ​​for individual typical nodal diameters, lacking a systematic assessment of the stability evolution within the complete traveling wave range (full nodal diameter), and failing to fully grasp the safety boundaries of the last-stage blade under all possible coupling modes. Summary of the Invention

[0008] The technical problem this invention aims to solve is to overcome the shortcomings of traditional models that neglect the effects of non-equilibrium condensation phase change. To this end, this invention proposes a method for predicting turbine blade flutter. This method utilizes the high efficiency of the influence coefficient method (ICM) to predict aerodynamic damping across the entire blade diameter, thereby accurately reflecting the true aeroelastic characteristics of the turbine's last-stage wet steam two-phase flow condition. This provides precise theoretical support for improving the reliability design of turbines for long-term service.

[0009] To achieve the above-mentioned objectives, the present invention specifically adopts the following technical solution:

[0010] In a first aspect, the present invention provides a method for predicting turbine blade flutter, comprising the following steps:

[0011] S1. Obtain the geometric and aerodynamic parameters of the turbine's last stage blade to be predicted, and use 3D modeling software to establish a 3D solid model of the turbine's last stage blade based on the geometric and aerodynamic parameters.

[0012] S2. Import the three-dimensional solid model into the finite element analysis software, perform tetrahedral meshing on it, and obtain the natural modes of the last stage blade of the steam turbine through modal analysis using the finite element analysis software. Extract the first three natural frequencies and corresponding mode shapes from them, and export the finite element mesh node information of the surface of the last stage blade of the steam turbine.

[0013] S3. Select an excitation frequency band within the first six natural frequencies obtained from the modal analysis of the last stage blade of the steam turbine. Select several excitation loads with adjacent interval frequencies within this excitation frequency band and apply them sequentially to the finite element mesh. Calculate the displacement response amplitude under different excitation loads through harmonic response analysis. Take the excitation frequency corresponding to the maximum displacement response amplitude as the main excitation frequency and the mode shape whose natural frequency is closest to the main excitation frequency as the main mode shape.

[0014] S4. Import the 3D solid model into the mesh generation software, perform hexahedral mesh generation on the last stage blades of the steam turbine, and expand the single-flow-domain mesh of the last stage blades of the steam turbine into a nine-channel flow field model.

[0015] S5. A modified wet steam model containing the mass generation rate equation, nucleation rate calculation equation, and droplet growth rate model is pre-constructed, and the modified wet steam model is loaded into the computational fluid domain of the nine-channel flow field model. Based on the correction effect of the average droplet radius on the droplet surface tension, the droplet surface tension coefficient in the nucleation rate calculation equation is corrected.

[0016] S6. Based on the modified wet steam model, the influence coefficient method is used to complete the numerical calculation of the unsteady flow field, and the modal force data of each blade channel in the last vibration cycle are obtained and output.

[0017] S7. Perform Discrete Fourier Transform on the output modal force data to extract its first harmonic components, thereby calculating the aerodynamic damping coefficient under the full pitch diameter. When the aerodynamic damping coefficient is greater than 0, it indicates that the flutter state of the last stage blade of the steam turbine is in the stable range. When the aerodynamic damping coefficient is less than 0, it indicates that the flutter state of the last stage blade of the steam turbine is in the unstable range, thus realizing the flutter stability prediction of the last stage moving blade of the steam turbine.

[0018] Based on the above scheme, each step can be implemented in the following preferred manner.

[0019] As a preferred embodiment of the first aspect above, in step S1, the geometric parameters include the blade pitch circle diameter, blade height, number of blades, and blade tip clearance; the aerodynamic parameters include the unit flow rate, total fluid inlet pressure, total fluid inlet temperature, static pressure at the fluid outlet, liquid phase volume fraction at the fluid inlet, and liquid phase diameter at the fluid inlet.

[0020] As a preferred embodiment of the first aspect mentioned above, in step S4, when performing hexahedral meshing on the three-dimensional solid model, the blade region of the last stage blade of the steam turbine adopts a HOH type mesh topology, the fluid inlet and fluid outlet of the last stage blade of the steam turbine both adopt H type orthogonal meshes, and the interface region between the fluid inlet and fluid outlet adopts an O type orthogonal mesh; local mesh refinement is performed in the blade tip clearance region, the blade leading edge and trailing edge, and the endwall boundary layer region.

[0021] As a preferred embodiment of the first aspect above, in step S5, the mass generation rate equation is expressed as:

[0022] ;

[0023]

[0024]

[0025] In the formula, For the rate of quality generation, This indicates the density of the liquid phase. Indicates the nucleation rate, The number of droplets per unit mass. Indicates the average droplet radius. This represents the partial derivative with respect to time. This represents the partial derivative with respect to the average droplet radius. Indicates the droplet growth rate. Indicates the critical droplet radius. This represents the corrected droplet surface tension coefficient. Represents the gas constant. Indicates the absolute temperature of steam. Indicates the supersaturation ratio. Indicates vapor pressure. This represents the equilibrium saturation pressure.

[0026] As a preferred embodiment of the first aspect above, in step S5, the equation for calculating the nucleation rate is expressed as:

[0027] ;

[0028]

[0029] In the formula, Indicates the evaporation-condensation coefficient. Indicates the temperature of the gas phase. Indicates the temperature of the liquid phase. This represents the density of the gas phase. Indicates the molecular mass of a single droplet. Represents the Boltzmann constant. This represents the non-isothermal correction factor. Represents the natural constant. Indicates the latent heat of vaporization. This indicates the specific heat ratio of steam.

[0030] As a preferred embodiment of the first aspect mentioned above, in step S5, the droplet growth rate model is as follows:

[0031] ;

[0032]

[0033]

[0034]

[0035] In the formula, Indicates the adjustment factor. Represents the dimensionless Knudsen number. The mean free path of the molecules, This represents the correction factor. The thermal conductivity of steam, Indicates supercooling. This represents the saturation temperature under the current pressure. For steam Prandtl number, Indicates the growth coefficient. This indicates the specific heat capacity of steam at constant pressure.

[0036] As a preferred embodiment of the first aspect mentioned above, in step S5, the corrected droplet surface tension coefficient is calculated as follows:

[0037] ;

[0038] ;

[0039] in, This represents the droplet surface tension coefficient before correction. Represents a constant variable. This indicates the critical temperature.

[0040] As a preferred embodiment of the first aspect mentioned above, step S6, which uses the influence coefficient method to complete the numerical calculation of the unsteady flow field, specifically includes: in the computational fluid dynamics software CFX, the principal mode shape determined in S3 is mapped to the blade surface nodes of the nine-channel flow field model through a three-dimensional linear interpolation method, so that the reference blade in the nine-channel flow field model undergoes simple harmonic vibration according to the principal mode shape and its corresponding natural angular frequency, while the other blades remain stationary, and unsteady numerical solutions are performed to obtain the modal force data of each blade channel within one vibration cycle.

[0041] In a second aspect, the present invention provides a computer-readable storage medium storing a computer program that, when executed by a processor, implements the turbine blade flutter prediction method as described in any of the solutions of the first aspect above.

[0042] Thirdly, the present invention provides a computer electronic device, which includes a memory and a processor;

[0043] The memory is used to store computer programs;

[0044] The processor is configured to, when executing the computer program, implement the turbine blade flutter prediction method as described in any of the first aspects above.

[0045] Compared with the prior art, the present invention has the following advantages:

[0046] The key point of the method provided by this invention lies in coupling the real physics of non-equilibrium wet steam with the influence coefficient method (ICM) for predicting flutter of long blades in the last stage of a steam turbine. This invention utilizes a developed adaptive non-equilibrium wet steam real physics model (i.e., a modified wet steam model) to replace the traditional single-phase equivalent model, fully considering the influence of phase change hysteresis, latent heat release, and condensation shock waves on the local pressure phase and amplitude, thus correcting the systematic deviation of aerodynamic work / damping at its source. Based on this, this invention employs structure-fluid modal consistency mapping (applying the dominant mode shape and natural angular frequency to the fluid wall with high fidelity) to ensure the aerodynamic-structure phase relationship remains undistorted, and combines this with end-period DFT extraction of the main harmonics. Furthermore, this invention leverages the finite-channel truncation and linear superposition of ICM to reconstruct the aerodynamic damping coefficients of the entire pitch diameter in a single multi-channel unsteady calculation, reducing the complexity of traditional TWM's successive recalculation by ND from "multiple full-cycle / multi-channel" to "single multi-channel," significantly reducing computational cost and time while avoiding incomplete traveling wave coverage. Meanwhile, the process naturally supports sensitive scanning of multiple parameters such as humidity, droplet size / number density, and back pressure, accurately outputting the negative damping risk zone, thereby providing highly reliable, full-coverage, and engineering-usable flutter prediction capabilities for the geometry optimization, end-zone reinforcement, and operation limit setting of the last-stage long blades. Attached Figure Description

[0047] Figure 1 This is a flowchart of the method of the present invention;

[0048] Figure 2 This is a schematic diagram illustrating the relationship between input frequency and displacement amplitude in this embodiment;

[0049] Figure 3 This is a schematic diagram of the structured mesh generation of the fluid domain in this embodiment;

[0050] Figure 4 This is a schematic diagram illustrating the blade numbering definition in this embodiment;

[0051] Figure 5 This is a schematic diagram of the periodically changing nine-channel modal forces in this embodiment;

[0052] Figure 6 This is a diagram showing the distribution of aerodynamic damping coefficients in this embodiment. Detailed Implementation

[0053] To make the above-mentioned objects, features, and advantages of the present invention more apparent and understandable, the specific embodiments of the present invention will be described in detail below with reference to the accompanying drawings. Many specific details are set forth in the following description to provide a thorough understanding of the present invention. However, the present invention can be practiced in many other ways different from those described herein, and those skilled in the art can make similar modifications without departing from the spirit of the present invention. Therefore, the present invention is not limited to the specific embodiments disclosed below. Technical features in the various embodiments of the present invention can be combined accordingly without mutual conflict.

[0054] In the description of this invention, it should be understood that the terms "first" and "second" are used only for descriptive purposes and should not be construed as indicating or implying relative importance or implicitly specifying the number of indicated technical features. Therefore, a feature defined with "first" and "second" may explicitly or implicitly include at least one of those features.

[0055] Traditional methods for predicting flutter in the last stage of steam turbines still use ideal gas or other single-phase flow with equivalent thermodynamic parameters as the working fluid and employ the Traveling Wave Mode (TWM) algorithm to predict aerodynamic damping under typical nodal diameters. However, long last-stage blades in steam turbines are more susceptible to aerodynamic excitation-induced flutter under wet steam conditions, and non-equilibrium phase transitions significantly exacerbate the uncertainty of their aerodynamic stability. Furthermore, the minimum number of flow channels required for modeling varies for different nodal diameters using the TWM method, sometimes even requiring modeling the entire circumference, and each calculation only yields results for one nodal diameter. Currently, based on Integrated Circuit Modeling (ICM), significant time costs can be reduced; that is, by performing unsteady calculations on a multi-blade channel, the aerodynamic damping coefficients for all nodal diameters can be obtained. Based on this, this invention designs a method for predicting flutter in steam turbine blades.

[0056] like Figure 1 As shown, in a preferred embodiment of the present invention, the turbine blade flutter prediction method includes the following steps S1 to S7. The specific implementation process of each step will be described in detail below.

[0057] S1. Obtain the geometric and aerodynamic parameters of the turbine's last stage blade to be predicted, and use 3D modeling software to establish a 3D solid model of the turbine's last stage blade based on the geometric and aerodynamic parameters.

[0058] It should be noted that in step S1 of the present invention, the geometric parameters include the blade pitch circle diameter, blade height, number of blades, and blade tip clearance; the aerodynamic parameters include the unit flow rate, total fluid inlet pressure, total fluid inlet temperature, static pressure at the fluid outlet, liquid phase volume fraction at the fluid inlet, and liquid phase diameter at the fluid inlet.

[0059] In step S1 of this embodiment, the mechanical damping settings such as the shroud are simplified using SolidWorks 3D modeling software. Basic geometric and aerodynamic parameters are then used to establish a 3D solid model of the turbine's last-stage moving blades to be analyzed. Specific design parameters are shown in Table 1.

[0060] Table 1. Basic Design Parameters of the Last Stage Moving Blades of Steam Turbines

[0061]

[0062] S2. Import the 3D solid model into the finite element analysis software, perform tetrahedral meshing on it, and obtain the natural modes of the last stage blade of the steam turbine through modal analysis using the finite element analysis software. Extract the first three natural frequencies and corresponding mode shapes from them, and export the finite element mesh node information of the surface of the last stage blade of the steam turbine.

[0063] In step S2 of this embodiment, the established three-dimensional solid model is imported into the finite element analysis software ANSYS. Its built-in meshing tool is used to generate a tetrahedral mesh, and material mechanical properties are assigned to the finite element mesh model, defining its material and element properties, as detailed in Table 2. Fixed support boundary constraints are applied to the blade root of the finite element mesh model, while the blade tip is free. Simultaneously, a centrifugal force load corresponding to the target rotational speed is applied. The natural modes of the turbine's last-stage blades are obtained using dynamic modal analysis in ANSYS. Dynamic modal analysis is performed on the finite element mesh model with applied constraints and centrifugal force load, extracting the first three frequencies and mode shapes. The finite element mesh node information of the turbine's last-stage blade surface is then exported, constructing the finite element mesh model.

[0064] In step S2 of this embodiment, the total number of finite element meshes is 84993; the centrifugal force load is set to 3000 r / min, which is closer to the actual working state of the last stage blade of the steam turbine.

[0065] Table 2. Parameters of the finite element model of the last stage moving blade structural domain of the steam turbine

[0066]

[0067] S3. Select an excitation frequency band within the first six natural frequencies obtained from the modal analysis of the last stage blade of the steam turbine. Select several excitation loads with adjacent interval frequencies within this excitation frequency band and apply them sequentially to the finite element mesh. Calculate the displacement response amplitude under different excitation loads through harmonic response analysis. Take the excitation frequency corresponding to the maximum displacement response amplitude as the main excitation frequency and the mode shape whose natural frequency is closest to the main excitation frequency as the main mode shape.

[0068] In step S3 of this embodiment, the first six natural frequencies of the turbine's last-stage blades are calculated based on modal analysis. The minimum frequency of the excitation band is set to 0Hz, and the maximum frequency is set to 1200Hz. Since the minimum frequency difference of the excitation band is the difference between the second and third natural frequencies, i.e., 130Hz, it is divided into 12 parts, with each 100Hz serving as an excitation load. The purpose is to generate a corresponding external load excitation response near each natural frequency. The relationship between the input frequency and the displacement amplitude is obtained sequentially based on the applied excitation loads, as follows: Figure 2 As shown, it can be seen that the frequencies corresponding to the maximum amplitude peaks in the excitation frequency responses in the X, Y, and Z directions are all close to the first-order mode frequencies. This indicates that during the vibration process of the last-stage turbine blade, the first-order bending mode is the main mode of the blade. Therefore, subsequent flutter studies will only focus on the first-order bending mode of the blade.

[0069] S4. Import the 3D solid model into the mesh generation software, perform hexahedral mesh generation on the last stage blades of the steam turbine, and expand the single-flow-domain mesh of the last stage blades of the steam turbine into a nine-channel flow field model.

[0070] It should be noted that in step S4 of this invention, when performing hexahedral meshing on the three-dimensional solid model, the blade region of the last stage blade of the steam turbine adopts a HOH type mesh topology, the fluid inlet and fluid outlet of the last stage blade of the steam turbine both adopt H type orthogonal meshes, and the interface region between the fluid inlet and fluid outlet adopts an O type orthogonal mesh; local mesh refinement is performed in the blade tip clearance region, the blade leading edge and trailing edge, and the endwall boundary layer region.

[0071] In step S4 of this embodiment, the commercial finite element software TurboGrid is used as the meshing software. TurboGrid uses a "HOH" type mesh topology to perform hexahedral meshing on the three-dimensional solid model. Both the fluid inlet and outlet use "H" type orthogonal meshes, while the interface region between the fluid inlet and outlet uses an "O" type orthogonal mesh to reduce mesh distortion at the leading edge, trailing edge, and tip clearance. Figure 3 As shown. In addition, due to various disturbances on the static wall surface, the hexahedral mesh is locally refined at the tip gap, leading and trailing edges, and endwall boundary layer to capture more complex flow phenomena.

[0072] Furthermore, in order to meet For the turbulence model requirements, in this embodiment, the y+ values ​​of the single-channel fluid domain mesh are all less than 30, the minimum angle is greater than 20°, and the mesh quality is greater than 0.5. The number of fluid domain meshes is 294,464. The single-channel mesh is then expanded to nine channels using ICEM software, resulting in a total of 2,985,180 fluid domain meshes.

[0073] S5. A modified wet steam model containing the mass generation rate equation, nucleation rate calculation equation, and droplet growth rate model is pre-constructed, and the modified wet steam model is loaded into the computational fluid domain of the nine-channel flow field model. Based on the correction effect of the average droplet radius on the droplet surface tension, the droplet surface tension coefficient in the nucleation rate calculation equation is corrected.

[0074] It should be noted that in step S5 of this invention, the mass generation rate equation is expressed as:

[0075] ;

[0076]

[0077]

[0078] In the formula, For the rate of quality generation, This indicates the density of the liquid phase. Indicates the nucleation rate, The number of droplets per unit mass. Indicates the average droplet radius. This represents the partial derivative with respect to time. This represents the partial derivative with respect to the average droplet radius. Indicates the droplet growth rate. Indicates the critical droplet radius. This represents the corrected droplet surface tension coefficient. Represents the gas constant. Indicates the absolute temperature of steam. Indicates the supersaturation ratio. Indicates vapor pressure. This represents the equilibrium saturation pressure.

[0079] It should be noted that in step S5 of this invention, the nucleation rate calculation equation is expressed as:

[0080] ;

[0081]

[0082] In the formula, Indicates the evaporation-condensation coefficient. Indicates the temperature of the gas phase. Indicates the temperature of the liquid phase. This represents the density of the gas phase. Indicates the molecular mass of a single droplet. Represents the Boltzmann constant. This represents the non-isothermal correction factor. Represents the natural constant. Indicates the latent heat of vaporization. This indicates the specific heat ratio of steam.

[0083] It should be noted that, in step S5 of this invention, the droplet growth rate model is as follows:

[0084] ;

[0085]

[0086]

[0087]

[0088] In the formula, Indicates the adjustment factor. Represents the dimensionless Knudsen number. The mean free path of the molecules, This represents the correction factor. The thermal conductivity of steam, Indicates supercooling. This represents the saturation temperature under the current pressure. For steam Prandtl number, Indicates the growth coefficient. This indicates the specific heat capacity of steam at constant pressure.

[0089] It should be noted that in step S5 of this invention, the corrected droplet surface tension coefficient is calculated as follows:

[0090] ;

[0091] ;

[0092] in, This represents the droplet surface tension coefficient before correction. Represents a constant variable. This indicates the critical temperature.

[0093] In step S5 of this embodiment, a modified wet steam model is constructed using an F source file written in FORTRAN and loaded into the computational fluid domain of the nine-channel flow field model. In the construction of the modified wet steam model, droplet growth is characterized by the average representative radius (i.e., the aforementioned average droplet radius), and the droplets are assumed to be ideally spherical and in an infinitely extending steam environment, thus ignoring boundary effects. Compared to the latent heat released during condensation, the heat capacity of small droplets is considered negligible. Under non-equilibrium condensation conditions, the mass generation rate in classical nucleation theory is defined as the sum of the mass increase caused by the formation of the critical droplet radius and the mass change generated during subsequent growth or decay of the droplet, used to comprehensively describe the contribution of the nucleation and growth processes to the overall mass change. The formula for calculating the mass generation rate is as described above and will not be repeated. Specifically, when the critical droplet radius is exceeded, the droplet will grow; when it is below this radius, the droplet will evaporate.

[0094] In nonequilibrium condensation studies, determining the nucleation rate and droplet growth rate are key factors affecting the accuracy of condensation characteristic predictions. For the nucleation rate, the modified wet steam model employs a modified Kantrowitz model. This model considers the influence of the temperature difference between steam and droplets on the phase change process, and the calculation results more accurately reflect the thermodynamic nonequilibrium characteristics of the actual condensation process, thus improving accuracy. The calculation method for the nucleation rate is as described above and will not be repeated. In this embodiment, the evaporation-condensation coefficient is set to 1, and the non-isothermal correction coefficient is used to introduce the temperature difference effect, thereby correcting the energy transfer process of the classical isothermal homogeneous nucleation theory.

[0095] In the exponential term of the above nucleation rate calculation equation, the corrected droplet surface tension coefficient is extremely sensitive to the nucleation rate; even a small change in it can lead to significant fluctuations in the condensation rate. To reasonably correct the nucleation rate calculation equation, this invention corrects the droplet surface tension coefficient in the equation based on the correction effect of the average droplet radius on the droplet surface tension. The calculation method for the corrected droplet surface tension coefficient is as described above and will not be repeated here.

[0096] The droplet surface tension coefficient before the correction was obtained from the relationship proposed by IAPWS-IF97. In this embodiment, the constant variable is taken as... The critical temperature is Experiments have verified that the corrected method for calculating the droplet surface tension coefficient is effective across the entire vapor-liquid saturation line starting from the triple point temperature.

[0097] In the field of nonequilibrium condensation research, accurate modeling of droplet growth rate is one of the core problems that urgently needs to be solved. The droplet growth process essentially involves two types of transport mechanisms: the first is the mass migration process from vapor to the droplet, and the second is the release of latent heat from the droplet to the surrounding vapor. To address the problem of large calculation deviations in droplet growth rate in traditional wet steam models, this embodiment designs a droplet growth rate model, the specific form of which is as described above and will not be repeated here. This model, while considering multiple effects of momentum and energy transfer, introduces Knudsen correction and non-isothermal correction, enabling it to more accurately reflect the droplet dynamics under nonequilibrium condensation conditions in numerical simulations.

[0098] In the above droplet growth rate model, and All of these are fitted parameters. It is the growth coefficient that is related to the condensation coefficient and the evaporation coefficient; This is an adjustment coefficient used to adjust the distance from the droplet to where the continuous process (relative to the free molecular process) occurs. The value is optimized by comparing with experimental data.

[0099] In this embodiment, The value is 9. The value is set to 1. The droplet growth rate model is embedded into the wet steam control equation, and the droplet generation rate in the nucleation theory is synchronously correlated to form a complete modified wet steam model. This provides reliable physical model support for subsequent numerical calculations of unsteady flow fields based on the influence coefficient method, and effectively improves the prediction accuracy of key parameters of the wet steam flow field.

[0100] S6. Based on the modified wet steam model, the influence coefficient method is used to complete the numerical calculation of the unsteady flow field, and the modal force data of each blade channel in the last vibration cycle are obtained and output.

[0101] It should be noted that, in step S6 of this invention, the numerical calculation process of the unsteady flow field using the influence coefficient method specifically includes: in the computational fluid dynamics software CFX, the principal mode shape determined in S3 is mapped to the blade surface nodes of the nine-channel flow field model using a three-dimensional linear interpolation method, so that the reference blade (i.e., blade 0) in the nine-channel flow field model is aligned with the principal mode shape and its corresponding natural angular frequency. Simple harmonic motion is performed while the other blades remain stationary. Unsteady numerical solutions are then obtained to acquire modal force data for each blade channel within one vibration cycle.

[0102] In step S6 of this embodiment, based on the above process, a modified wet steam model is compiled and loaded, and boundary conditions are loaded into the nine-channel flow field model established in step S4. The inlet and outlet refer to Table 1, and a cyclic symmetry plane is loaded on the two outer channels. Motion mesh boundary conditions are given on the blade surface. The working fluid uses steam3v from the IAPWS-IF97 library of CFX software as the gas phase and steam3l as the liquid phase, and is compared with an ideal gas using equivalent thermodynamic properties. In this embodiment, the specific heat ratio of the ideal gas with equivalent thermodynamic properties is... Specific heat capacity of steam at constant pressure The dynamic viscosity is Thermal conductivity is Using a steady field without a moving mesh as the initial condition, unsteady calculations were performed on two different models using the influence coefficient method.

[0103] In the technical framework for solving unsteady aerodynamic forces, the influence coefficient method is based on two core assumptions: the first is the finite propagation assumption, which limits the circumferential propagation of disturbances caused by blade vibration, meaning that vibration disturbances can only affect local blades on both sides of the circumference, rather than propagating unbounded along the circumference; the second is the linear superposition assumption, which states that unsteady pressures at different phase angles can be obtained by linearly superimposing the influence coefficients. When using the influence coefficient method to solve unsteady aerodynamic forces, it is necessary to specify a single blade to vibrate at a specific frequency and mode shape, extract aerodynamic data from different blade surfaces, and then perform linear superposition calculations using the following formula to obtain the total unsteady pressure of different traveling waves. :

[0104] ;

[0105] in, For the phase difference of the blades, Indicates the blade number. Represents the imaginary unit. Indicates the number of channels used. Symbols indicate direction. This indicates the unsteady pressure caused by blade vibration. This represents the unsteady pressure of a single blade, with the subscripts TW and IC corresponding to the traveling wave and influence coefficient, respectively.

[0106] No. Modal forces of each blade Defined as:

[0107] ;

[0108] in, For the first Instantaneous static pressure on each blade Indicates the first One leaf, , , These are the mode shape vectors in the x, y, and z directions, respectively. , , These represent the unit normals of the blade surface in the x, y, and z directions, respectively. Indicates the surface of the blade. This represents the integral over the surface area of ​​the blade.

[0109] In this embodiment, as Figure 4 As shown, the blade numbering is defined as follows: the vibrating blade is the reference blade, with a number of 0. Blades upstream of the reference blade in the direction of rotation are numbered starting from 1, gradually increasing in number as they move away from the reference blade. Blades downstream in the direction of rotation are numbered starting from -1, gradually decreasing in number as they move away from the reference blade. For the nine-channel computational domain, the blade number ranges from -4 to +4.

[0110] In the nine-channel computational domain of this embodiment, only the middle blade (0) undergoes simple harmonic motion according to its first-order modal natural frequency and mode shape, with a specified small amplitude of 0.001 meters, approximately 1% of the blade tip chord length. The remaining blades remain stationary, and their motion period is the reciprocal of the first-order modal frequency. Each period is divided into 100 time steps, and a total of 8 periods are calculated. At each time step, a modified wet steam model is applied to solve the Navier-Stokes equations. The output is set to output the modal force and corresponding nodal displacement information of each blade channel at each time step. The final result is as follows: Figure 5 The nine-channel modal forces are shown to change periodically.

[0111] S7. Perform Discrete Fourier Transform on the output modal force data to extract its first harmonic components, thereby calculating the aerodynamic damping coefficient under the full pitch diameter. When the aerodynamic damping coefficient is greater than 0, it indicates that the flutter state of the last stage blade of the steam turbine is in the stable range. When the aerodynamic damping coefficient is less than 0, it indicates that the flutter state of the last stage blade of the steam turbine is in the unstable range, thus realizing the flutter stability prediction of the last stage moving blade of the steam turbine.

[0112] In step S7 of this embodiment, the modal force data is subjected to discrete Fourier transform to obtain the first harmonic components of the modal force:

[0113] ;

[0114]

[0115] in, This indicates that the vibration of the reference blade affects the first... Modal forces generated by each blade The total number of discrete time steps within one oscillation cycle. Set to 100, For discrete time steps, Let a discrete time step represent the time. Indicates the first The blades vibrate due to the vibration of the reference blade. Modal forces generated at any time A vibration period can be represented by the reciprocal of the frequency in the first-order mode shape.

[0116] In this embodiment, further definition Let be the tip phase angle between leaf 1 and leaf 0 (reference leaf), then the... The phase difference between blade 0 and blade 1 is Therefore, refer to the blade because of the first The modal force generated by the vibration of each blade can be expressed as: Assume the blade has a modal displacement of 0. Vibration, then from the first The work done by the modal forces caused by the vibration of a blade during one vibration cycle on blade 0 for:

[0117]

[0118]

[0119] in, Indicates modal displacement, Indicates taking the real part, Indicates the number of leaves. This indicates the node diameter. Here, the number of blades is set to 65.

[0120] By summing the work done by all the blades, we can obtain the total aerodynamic work of the traveling wave vibration corresponding to the phase difference. :

[0121] ;

[0122] This leads to the aerodynamic damping coefficient. The calculation formula:

[0123] ;

[0124] in, This represents the natural angular frequency.

[0125] Therefore, the aerodynamic damping coefficient under full-diameter conditions can be obtained through the above calculations. Based on this, further analysis of the aerodynamic damping coefficient under full-diameter conditions can be conducted to establish the correspondence between aerodynamic damping and diameter and wet steam operating conditions, such as... Figure 6 As shown, this method enables the prediction of flutter stability and identification of critical pitch diameters for the last-stage moving blades of a steam turbine. This approach aims to overcome the shortcomings of traditional models that neglect non-equilibrium condensation phase change effects, achieving aerodynamic damping prediction across the entire pitch diameter. This accurately reflects the true aeroelastic characteristics of the last-stage moving blades under wet steam two-phase flow conditions, significantly reducing computational time costs and improving the physical reliability and reproducibility of the predictions. This provides precise theoretical support for improving the reliability design of steam turbines for long-term service.

[0126] It is understood that the turbine blade flutter prediction method described in S1 to S7 above can essentially be implemented by a computer program. Therefore, based on the same inventive concept, another preferred embodiment of the present invention also provides a computer electronic device corresponding to the turbine blade flutter prediction method provided in the above embodiments, which includes a memory and a processor;

[0127] The memory is used to store computer programs;

[0128] The processor is used to implement the turbine blade flutter prediction method in the above embodiments when executing the computer program.

[0129] Furthermore, the logical instructions in the aforementioned memory can be implemented as software functional units and, when sold or used as independent products, can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of the present invention, in essence, or the part that contributes to the prior art, or a portion of the technical solution, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute all or part of the steps of the methods described in the various embodiments of the present invention.

[0130] Therefore, based on the same inventive concept, another preferred embodiment of the present invention also provides a computer-readable storage medium corresponding to the turbine blade flutter prediction method provided in the above embodiments. The storage medium stores a computer program, which, when executed by a processor, can realize the turbine blade flutter prediction method in the above embodiments.

[0131] It is understood that the aforementioned storage media may include random access memory (RAM) or non-volatile memory (NVM), such as at least one disk storage device. Furthermore, the storage media may also be various media capable of storing program code, such as USB flash drives, external hard drives, magnetic disks, or optical discs.

[0132] It is understood that the processors mentioned above can be general-purpose processors, including central processing units (CPUs), network processors (NPs), etc.; they can also be digital signal processors (DSPs), application-specific integrated circuits (ASICs), field-programmable gate arrays (FPGAs), or other programmable logic devices, discrete gate or transistor logic devices, or discrete hardware components.

[0133] The embodiments described above are merely preferred embodiments of the present invention and are not intended to limit the invention. Those skilled in the art can make various changes and modifications without departing from the spirit and scope of the invention. Therefore, all technical solutions obtained through equivalent substitution or transformation fall within the protection scope of the present invention.

Claims

1. A method of flutter prediction for a turbine blade, characterized by, Includes the following steps: S1. Obtain the geometric and aerodynamic parameters of the turbine's last stage blade to be predicted, and use 3D modeling software to establish a 3D solid model of the turbine's last stage blade based on the geometric and aerodynamic parameters. S2. Import the three-dimensional solid model into the finite element analysis software, perform tetrahedral meshing on it, and obtain the natural modes of the last stage blade of the steam turbine through modal analysis using the finite element analysis software. Extract the first three natural frequencies and corresponding mode shapes from them, and export the finite element mesh node information of the surface of the last stage blade of the steam turbine. S3. Select an excitation frequency band within the first six natural frequencies obtained from the modal analysis of the last stage blade of the steam turbine. Select several excitation loads with adjacent interval frequencies within this excitation frequency band and apply them sequentially to the finite element mesh. Calculate the displacement response amplitude under different excitation loads through harmonic response analysis. Take the excitation frequency corresponding to the maximum displacement response amplitude as the main excitation frequency and the mode shape whose natural frequency is closest to the main excitation frequency as the main mode shape. S4. Import the 3D solid model into the mesh generation software, perform hexahedral mesh generation on the last stage blades of the steam turbine, and expand the single-flow-domain mesh of the last stage blades of the steam turbine into a nine-channel flow field model. S5. A modified wet steam model containing the mass generation rate equation, nucleation rate calculation equation, and droplet growth rate model is pre-constructed, and the modified wet steam model is loaded into the computational fluid domain of the nine-channel flow field model. Based on the correction effect of the average droplet radius on the droplet surface tension, the droplet surface tension coefficient in the nucleation rate calculation equation is corrected. S6. Based on the modified wet steam model, the influence coefficient method is used to complete the numerical calculation of the unsteady flow field, and the modal force data of each blade channel in the last vibration cycle are obtained and output. S7. Perform Discrete Fourier Transform on the output modal force data to extract its first harmonic components, thereby calculating the aerodynamic damping coefficient under the full pitch diameter. When the aerodynamic damping coefficient is greater than 0, it indicates that the flutter state of the last stage blade of the steam turbine is in the stable range. When the aerodynamic damping coefficient is less than 0, it indicates that the flutter state of the last stage blade of the steam turbine is in the unstable range, thus realizing the flutter stability prediction of the last stage moving blade of the steam turbine.

2. The turbine blade flutter prediction method of claim 1, wherein, In step S1, the geometric parameters include blade pitch circle diameter, blade height, number of blades, and blade tip clearance; the aerodynamic parameters include unit flow rate, total fluid inlet pressure, total fluid inlet temperature, static pressure at the fluid outlet, liquid phase volume fraction at the fluid inlet, and liquid phase diameter at the fluid inlet.

3. The turbine blade flutter prediction method of claim 1, wherein, In step S4, when performing hexahedral meshing on the three-dimensional solid model, the blade region of the last stage turbine blade adopts the HOH type mesh topology, the fluid inlet and fluid outlet of the last stage turbine blade both adopt the H type orthogonal mesh, and the interface region between the fluid inlet and fluid outlet adopts the O type orthogonal mesh; local mesh refinement is performed in the blade tip clearance region, the blade leading edge and trailing edge, and the endwall boundary layer region.

4. The turbine blade flutter prediction method of claim 1, wherein, In step S5, the mass generation rate equation is expressed as: ； ； ； In the formula, For the rate of quality generation, This indicates the density of the liquid phase. Indicates the nucleation rate, The number of droplets per unit mass. Indicates the average droplet radius. This represents the partial derivative with respect to time. This represents the partial derivative with respect to the average droplet radius. Indicates the droplet growth rate. Indicates the critical droplet radius. This represents the corrected droplet surface tension coefficient. Represents the gas constant. Indicates the absolute temperature of steam. Indicates the supersaturation ratio. Indicates vapor pressure. This represents the equilibrium saturation pressure.

5. The turbine blade flutter prediction method as described in claim 4, characterized in that, In step S5, the equation for calculating the nucleation rate is expressed as: ； ； In the formula, Indicates the evaporation-condensation coefficient. Indicates the temperature of the gas phase. Indicates the temperature of the liquid phase. This represents the density of the gas phase. Indicates the molecular mass of a single droplet. Represents Boltzmann's constant. This represents the non-isothermal correction factor. Represents the natural constant. Indicates the latent heat of vaporization. This indicates the specific heat ratio of steam.

6. The turbine blade flutter prediction method as described in claim 5, characterized in that, In step S5, the droplet growth rate model is as follows: ； ； ； ； In the formula, Indicates the adjustment factor. Represents the dimensionless Knudsen number. The mean free path of the molecules, This represents the correction factor. The thermal conductivity of steam, Indicates supercooling. This represents the saturation temperature under the current pressure. For steam Prandtl number, Indicates the growth coefficient. This indicates the specific heat capacity of steam at constant pressure.

7. The turbine blade flutter prediction method as described in claim 6, characterized in that, In step S5, the corrected droplet surface tension coefficient is calculated as follows: ； ； in, This represents the droplet surface tension coefficient before correction. Represents a constant variable. This indicates the critical temperature.

8. The turbine blade flutter prediction method as described in claim 1, characterized in that, In step S6, the numerical calculation process of the unsteady flow field using the influence coefficient method specifically includes: in the computational fluid dynamics software CFX, the principal mode shape determined in S3 is mapped to the blade surface nodes of the nine-channel flow field model through a three-dimensional linear interpolation method, so that the reference blade in the nine-channel flow field model undergoes simple harmonic vibration according to the principal mode shape and its corresponding natural angular frequency, while the other blades remain stationary, and unsteady numerical solutions are performed to obtain the modal force data of each blade channel in one vibration cycle.

9. A computer-readable storage medium, characterized in that, The storage medium stores a computer program, which, when executed by a processor, implements the turbine blade flutter prediction method as described in any one of claims 1 to 8.

10. A computer electronic device, characterized in that, Including memory and processor; The memory is used to store computer programs; The processor is configured to implement the turbine blade flutter prediction method as described in any one of claims 1 to 8 when executing the computer program.

Images