Nonlinear Blade Flutter Evaluation Method Based on Accumulated Work Measurement and Fatigue Life

Abstract

The present invention relates to a nonlinear blade flutter assessment method based on accumulated work measurement and fatigue life. The method obtains the maximum allowable blade vibration amplitude based on fatigue life as a target. Taking into account the blade nonlinearity caused by excessive amplitude in simple harmonic vibration, the blade is vibrated at the maximum allowable vibration amplitude. The harmonic balance method is then used to solve the unsteady flow, obtaining aerodynamic accumulated work. The positive or negative value of the accumulated work is then used to determine whether the blade flutter characteristics meet the design safety requirements. This method uses accumulated work measurement and fatigue life to determine whether blade vibrations are at risk of exceeding amplitude limits, effectively resolving the issue of conservatism in flutter linear analysis and improving the accuracy of blade flutter assessment.

Description

Technical Field

[0001] The present invention relates to the technical field of impeller machinery design, and in particular to a method for evaluating blade flutter by using accumulated work measurement. Background Art

[0002] Turbines are essential components of aircraft engine cores, and blades are crucial components. Blade flutter is a self-excited vibration, where the unsteady aerodynamic forces that cause blade oscillations originate from blade oscillations. It is a phenomenon known as aeroelastic instability. When flutter occurs, the blade vibration amplitude initially increases exponentially. Without adequate mechanical damping, the blade vibration can quickly exceed the limit and lead to blade fracture. However, with adequate mechanical damping, the blade vibration may enter a limit cycle state with a constant amplitude, never exceeding the maximum allowable vibration amplitude.

[0003] After flutter occurs, as the blade vibration amplitude increases, the blade vibration enters nonlinear fluid-solid coupling vibration and may enter limit cycle vibration. However, as long as the vibration amplitude does not exceed the limit, this flutter does not cause any substantial harm.

[0004] At present, the general flutter analysis method is a linear analysis method, which can only determine whether flutter occurs. That is, this flutter analysis method is relatively conservative. It does not focus on whether the vibration amplitude will increase to an excessive limit after flutter occurs, and thus cannot determine, let alone timely determine whether the blade design meets safety requirements. Summary of the Invention

[0005] The purpose of the present invention is to avoid the shortcomings of the prior art and provide a method for obtaining the maximum allowable vibration amplitude of the blade based on the fatigue life as the goal, and in the numerical analysis of aerodynamic damping, the blade is made to vibrate according to the maximum allowable vibration amplitude, and then the unsteady flow is solved by the harmonic balance method, that is, the flow nonlinearity caused by the excessive vibration amplitude of the blade is taken into account, and the aerodynamic accumulated work is obtained. Finally, the positive or negative value of the accumulated work is used to judge whether the flutter characteristics of the blade meet the requirements, thereby establishing a nonlinear blade flutter evaluation method based on accumulated work measurement and fatigue life.

[0006] To achieve the above objectives, the present invention adopts a technical solution: a nonlinear blade flutter assessment method based on accumulated work measurement and fatigue life, comprising the following steps:

[0007] Step 1: Establish the fluid domain calculation grid of the blade, perform steady-state calculation of the blade flow field, obtain the characteristic lines of flow rate and pressure ratio, flow rate and efficiency at different blade speeds, and thus determine the calculation conditions of blade flutter; at the same time, based on structural statics and dynamics analysis, obtain the modal vibration shapes corresponding to at least the first three vibration modes of the blade vibration and the natural angular frequency ;

[0008] Step 2: Any first-order mode shape Interpolate to the calculation grid of the fluid domain of the blade, and calculate the maximum allowable vibration amplitude of the blade based on the maximum allowable alternating stress that meets the high-cycle fatigue life of the blade ;

[0009] Step 3: Under the calculation conditions, let the blade vibrate in any mode , natural angular frequency and maximum permissible vibration amplitude Perform simple harmonic oscillation;

[0010] Under the premise of considering the blade nonlinearity caused by the excessive amplitude in the simple harmonic vibration of the blade, the fundamental frequency and harmonic information of the unsteady flow caused by the simple harmonic vibration of the blade are specified, and the unsteady blade flutter calculation is initialized based on the steady calculation results of the blade flow field. Thus, the harmonic balance method is used to analyze the unsteady flow caused by blade vibration and the unsteady pressure on the blade surface is obtained. and blade vibration speed ;

[0011] Finally, the unsteady pressure on the blade surface during one vibration cycle is and blade vibration speed Integrate to get the accumulated work;

[0012] Step 4: Repeat steps 2 to 3 to calculate the accumulated work corresponding to each vibration mode and evaluate the blade flutter. Specifically:

[0013] If the accumulated work is positive, the amount of work done by the airflow on the blade increases during the cycle, and the amplitude will continue to increase. There is a risk of exceeding the amplitude limit, and the fatigue life requirements of the corresponding vibration mode will not be met;

[0014] If the accumulated work is negative, it will not exceed the maximum allowable vibration amplitude, so there is no risk of exceeding the amplitude limit, and the fatigue life requirements of the corresponding vibration mode are met.

[0015] Furthermore, in the steady-state calculation of the blade flow field, the FANS equation for determining the calculation conditions of the blade flutter is solved, and the FANS equation is expressed in a cylindrical coordinate system as:

[0016]

[0017] Where: is physical time, is the conserved variable within the unit volume of the flow field, The expression is:

[0018] ,

[0019] is the fluid density, 、 、 are the velocity components of the fluid in the axial, circumferential and radial directions, is the total energy per unit mass of fluid; are the axial, circumferential and radial coordinate positions respectively, are the velocity components of the fluid domain calculation grid in the axial, circumferential and radial directions, 、 、 They represent the axial, circumferential, and radial fluxes generated by the movement of the computational grid in the fluid domain;

[0020] are the convective fluxes of the fluid in the axial, circumferential and radial directions respectively. The convective flux represents the physical quantity transported through a unit area in unit time. The expression formula is:

[0021] ,

[0022] is the total enthalpy per unit mass of fluid, is the static pressure of the fluid;

[0023] are the viscous fluxes of the fluid in the axial, circumferential and radial directions respectively. The viscous flux is used to describe the diffusion transport of momentum and energy caused by viscous stress and heat conduction. The expression is:

[0024]

[0025]

[0026] ,

[0027] Where, is the normal stress between fluids, is the shear stress between fluids, is the heat flow between fluids, 、 、 The heat flow between the fluids is Components in the axial, circumferential, and radial directions;

[0028] is the source term, which is used to describe the influence of the internal or external flow system on the conserved quantity. The expression is:

[0029] ,

[0030] After solving the FANS equation, the characteristic lines of flow rate and pressure ratio, flow rate and efficiency at different blade speeds are obtained, and the calculation conditions are selected according to the characteristic lines for flutter analysis;

[0031] The calculation conditions include design conditions, congestion conditions and near-stall conditions.

[0032] Furthermore, in step 1, the modal vibration modes corresponding to at least the first three vibration modes of the blade vibration are obtained. and the natural angular frequency , specifically:

[0033] Structural dynamics equations based on blade vibration:

[0034] ,

[0035] Where: is the mass matrix, is the damping matrix, is the stiffness matrix, f is the aerodynamic force acting on the blade surface, is the vibration displacement of the blade, for The first derivative of is the vibration speed of the blade; is the acceleration of the blade;

[0036] In modal analysis, aerodynamic forces are not considered. f With the damping matrix C, we get formula 1:

[0037] ,

[0038] make Formula 2:

[0039] ,

[0040] Where: is the natural angular frequency of the blade, t is the physical time, and after substituting Formula 2 into Formula 1, we get Formula 3:

[0041] ,

[0042] Formula 3 is the characteristic equation. Solving Formula 3 will give the modal vibration shapes corresponding to each vibration mode of the blade. and the natural angular frequency , and the mode shapes The corresponding modal stress and modal displacement .

[0043] Furthermore, the step 2 is specifically as follows:

[0044] The maximum allowable alternating stress that satisfies the high cycle fatigue life of the blade is obtained from the Goodman diagram , and then, the maximum allowable vibration amplitude of the blade is obtained L for:

[0045] ,

[0046] Where, is the modal stress, is the modal displacement.

[0047] Furthermore, in step three:

[0048] Specify the fundamental frequency of the unsteady flow caused by blade vibration, which is the natural frequency of the blade :

[0049] ,

[0050] Where, is the natural angular frequency;

[0051] At the same time, specify the inter-leaf phase angle , interleaf phase angle That is, the phase difference between the blades when vibrating in the traveling wave mode, which can be expressed as:

[0052] ,

[0053] Where: is the vibration pitch diameter of the blade, is the number of leaves;

[0054] when When it is an even number, The value is [- / 2, / 2), where - / 2 and / 2 difference , only one of them is taken during calculation;

[0055] when When is an odd number, The value is [-( -1) / 2, ( -1) / 2];

[0056] In addition, to resolve the unsteady flow caused by blade vibration, the specified harmonic information is:

[0057] ,

[0058] in, is the harmonic number, and is a positive integer greater than 1.

[0059] Furthermore, the unsteady blade flutter calculation is initialized based on the steady calculation results of the blade flow field in step 3, and then the unsteady flow of the blade caused by the blade vibration is analyzed by the harmonic balance method to obtain the unsteady pressure on the blade surface corresponding to the unsteady flow of the blade. and blade vibration speed , wherein the harmonic balance method formula is:

[0060] ,

[0061] Where: is the variable of the unsteady flow field around the blade, is the time-averaged quantity of the unsteady flow field, is the physical time, n is the harmonic number used to analyze the unsteady flow field, i is the order corresponding to the different harmonics, and For the i The Fourier coefficients corresponding to the subharmonics are, It is i The angular frequency of the unsteady flow field corresponding to the subharmonic; the harmonic balance method is used to express the flow field variables of the unsteady flow as the time-averaged quantity plus the pulsating quantity, and approximate the pulsating quantity with the Fourier series;

[0062] Harmonic balance method is used to represent the unsteady pressure on the blade surface The expression formula is:

[0063] ,

[0064] Where: is the time-averaged unsteady pressure on the blade surface;

[0065] Similarly, the vibration speed of the blade is obtained The expression is:

[0066] ,

[0067] Where: is the time-averaged value of the blade vibration velocity.

[0068] Further, the accumulated work described in step three for:

[0069] ,

[0070] Where: is the vibration speed of the blade, is the fluid domain mesh normal vector on the blade surface, is the differential of the leaf surface area, is the blade vibration period, dt is the physical time differential, and S is the entire blade surface area.

[0071] Furthermore, the above analysis based on structural statics and dynamics obtains at least the first three modal vibration modes that are prone to flutter. and the natural angular frequency .

[0072] The beneficial effects of the present invention are as follows: the evaluation method proposed in the present invention obtains the maximum allowable vibration amplitude of the blade according to the fatigue life target, performs simple harmonic vibration with the maximum allowable amplitude, solves the unsteady flow using the harmonic balance method, and calculates the aerodynamic accumulated work, that is, the flow nonlinearity caused by the excessive blade vibration amplitude is taken into account. When the flow nonlinearity caused by the blade vibration is strong, high-order harmonics are introduced to deal with the unsteady flow, and then, based on the positive or negative accumulated work, it is judged whether the blade flutter characteristics meet the design safety requirements and whether there is a risk of the blade vibration exceeding the amplitude limit, which effectively solves the problem of over-conservative flutter linear analysis. BRIEF DESCRIPTION OF THE DRAWINGS

[0073] Figure 1 Schematic diagram of the method flow of the present invention;

[0074] Figure 2 is the calculation domain of the blade simulation in the specific example described in the present invention;

[0075] Figure 3 is a schematic diagram of the blade modal displacement in the specific calculation example of the present invention;

[0076] Figure 4 is a schematic diagram of blade modal stress in a specific example of the present invention;

[0077] Figure 5 Schematic diagram of the accumulated work density distribution on the blade surface when the blade vibrates in the first-order mode and the vibration pitch diameter is 0 in a specific example of the present invention. DETAILED DESCRIPTION

[0078] The principles and features of the present invention are described below with reference to the accompanying drawings. The examples given are only used to explain the present invention and are not used to limit the scope of the present invention.

[0079] In order to achieve the above object, the present invention provides the following specific implementation methods:

[0080] Example 1: Figure 1 As shown, a nonlinear blade flutter assessment method based on accumulated work measurement and fatigue life includes the following steps:

[0081] S01. Establish a computational grid for the fluid domain of the blade and perform steady-state calculations of the blade flow field. In the steady-state calculations of the blade flow field, solve the FANS equations used to determine the computational conditions of blade flutter and express the FANS equations in the cylindrical coordinate system as follows:

[0082]

[0083] Where: is physical time, is the conserved variable within the unit volume of the flow field, The expression is:

[0084] ,

[0085] is the fluid density, 、 、 are the velocity components of the fluid in the axial, circumferential and radial directions, is the total energy per unit mass of fluid; are the axial, circumferential and radial coordinate positions respectively, are the velocity components of the fluid domain calculation grid in the axial, circumferential and radial directions, 、 、 They represent the axial, circumferential, and radial fluxes generated by the movement of the computational grid in the fluid domain;

[0086] are the convective fluxes of the fluid in the axial, circumferential and radial directions respectively. The convective flux represents the physical quantity transported through a unit area in unit time. The expression formula is:

[0087] ,

[0088] is the total enthalpy per unit mass of fluid, is the static pressure of the fluid;

[0089] are the viscous fluxes of the fluid in the axial, circumferential and radial directions respectively. The viscous flux is used to describe the diffusion transport of momentum and energy caused by viscous stress and heat conduction. The expression is:

[0090]

[0091]

[0092] ,

[0093] Where, is the normal stress between fluids, is the shear stress between fluids, is the heat flow between fluids, 、 、 The heat flow between the fluids is Components in the axial, circumferential, and radial directions;

[0094] is the source term, which is used to describe the influence of the internal or external flow system on the conserved quantity. The expression is:

[0095] ,

[0096] After solving the FANS equation, the characteristic lines of flow rate and pressure ratio, and flow rate and efficiency at different blade speeds are obtained, thereby determining the calculation conditions for blade flutter. The calculation conditions can then be selected based on the characteristic lines for subsequent flutter analysis.

[0097] The calculation conditions include design conditions, congestion conditions and near-stall conditions.

[0098] S02. Based on structural statics and dynamics analysis, obtain the modal vibration shapes corresponding to at least the first three vibration modes of the blade vibration. and the natural angular frequency , specifically:

[0099] Structural dynamics equations based on blade vibration:

[0100] ,

[0101] Where: is the mass matrix, is the damping matrix, is the stiffness matrix, f is the aerodynamic force acting on the blade surface, is the vibration displacement of the blade, for The first derivative of is the vibration speed of the blade; is the acceleration of the blade;

[0102] In modal analysis, aerodynamic forces are not considered. f With the damping matrix C, we get formula 1:

[0103] ,

[0104] make Formula 2:

[0105] ,

[0106] Where: is the natural angular frequency of the blade, t is the physical time, and after substituting Formula 2 into Formula 1, we get Formula 3:

[0107] ,

[0108] Formula 3 is the characteristic equation. Solving Formula 3 will give the modal vibration shapes corresponding to each vibration mode of the blade. and the natural angular frequency , and the mode shapes The corresponding modal stress and modal displacement .

[0109] S03, any first-order mode vibration shape Interpolate to the computational grid of the blade's fluid domain and obtain the maximum allowable alternating stress that satisfies the blade's high-cycle fatigue life from the Goodman diagram , and then, the maximum allowable vibration amplitude of the blade is obtained L for:

[0110] ,

[0111] Where, is the modal stress, is the modal displacement.

[0112] S04. Under the calculation conditions, let the blade vibrate in any order mode. , natural angular frequency and maximum permissible vibration amplitude Perform simple harmonic vibration and, taking into account the blade nonlinearity caused by the large amplitude in the simple harmonic vibration of the blade, specify the fundamental frequency of the unsteady flow caused by the blade vibration, which is the natural frequency of the blade. :

[0113] ,

[0114] Where, is the natural angular frequency;

[0115] At the same time, specify the inter-leaf phase angle , interleaf phase angle That is, the phase difference between the blades when vibrating in the traveling wave mode, which can be expressed as:

[0116] ,

[0117] Where: is the vibration pitch diameter of the blade, is the number of leaves;

[0118] when When it is an even number, The value is [- / 2, / 2), where - / 2 and / 2 difference , only one of them is taken during calculation;

[0119] when When is an odd number, The value is [-( -1) / 2, ( -1) / 2];

[0120] In addition, to resolve the unsteady flow caused by blade vibration, the specified harmonic information is:

[0121] ,

[0122] in, is the harmonic number, and is a positive integer greater than 1.

[0123] S05. Initialize the unsteady blade flutter calculation based on the steady calculation results of the blade flow field, and then use the harmonic balance method to analyze the unsteady flow of the blade caused by the blade vibration to obtain the unsteady pressure on the blade surface corresponding to the unsteady flow of the blade. and blade vibration speed , wherein the harmonic balance method formula is:

[0124] ,

[0125] Where: is the variable of the unsteady flow field around the blade, is the time-averaged quantity of the unsteady flow field, is the physical time, n is the harmonic number used to analyze the unsteady flow field, i is the order corresponding to the different harmonics, and For the i The Fourier coefficients corresponding to the subharmonics are, It is i The angular frequency of the unsteady flow field corresponding to the subharmonic; the harmonic balance method is used to express the flow field variables of the unsteady flow as the time-averaged quantity plus the pulsating quantity, and approximate the pulsating quantity with the Fourier series;

[0126] Harmonic balance method is used to represent the unsteady pressure on the blade surface The expression formula is:

[0127] ,

[0128] Where: is the time-averaged unsteady pressure on the blade surface;

[0129] Similarly, the vibration speed of the blade is obtained The expression is:

[0130] ,

[0131] Where: is the time-averaged value of the blade vibration velocity.

[0132] S06. Unsteady pressure on the blade surface during one vibration cycle and blade vibration speed Integrate to get the accumulated work for:

[0133] ,

[0134] Where, is the vibration speed of the blade, is the fluid domain mesh normal vector on the blade surface, is the differential of the leaf surface area, is the blade vibration period, dt is the physical time differential, and S is the entire blade surface area.

[0135] S07. Repeat S03 to S06 to calculate the accumulated work corresponding to each vibration mode and evaluate the blade flutter. Specifically,

[0136] If the accumulated work is positive, the amount of work done by the airflow on the blade increases during the cycle, and the amplitude will continue to increase. There is a risk of exceeding the amplitude limit, and the fatigue life requirements of the corresponding vibration mode will not be met;

[0137] If the accumulated work is negative, it will not exceed the maximum allowable vibration amplitude, so there is no risk of exceeding the amplitude limit, and the fatigue life requirements of the corresponding vibration mode are met.

[0138] Example 2: Same as Example 1, except that in S02, based on the structural statics and dynamics analysis, at least the first three modal vibration modes that are prone to flutter are obtained. and the natural angular frequency .

[0139] In order to further illustrate the solution and beneficial effects of the present invention, a specific calculation example using NASA Rotor 67 is provided:

[0140] NASA Rotor 67 is a fan rotor in a multi-stage compressor designed and tested by NASA in the late 1980s. It has 22 blades per cycle and a design speed of 16043 rpm.

[0141] like Figure 2-4 As shown in FIG, a nonlinear blade flutter assessment method based on blade aerodynamic accumulated work measurement includes the following steps:

[0142] Step 1: Perform steady-state calculations on the blade flow field to obtain the flow-pressure ratio and flow-efficiency characteristic lines at different speeds, and determine the calculation conditions for blade flutter. Specifically:

[0143] Generate the computational mesh of the fluid domain of the blade, such as Figure 2 As shown. Flow field calculation requires solving the FANS equation, which is expressed in the cylindrical coordinate system as follows:

[0144]

[0145] Where: is a conserved variable, is the source term, is physical time, are the axial, circumferential and radial convection fluxes respectively, are the axial, circumferential, and radial viscous fluxes, respectively; are the velocity components of the fluid domain calculation grid in the axial, circumferential and radial directions, 、 、 They represent the axial, circumferential, and radial fluxes generated by the movement of the computational grid in the fluid domain;

[0146] After the flow field calculation was completed, the characteristic curve of the Rotor 67 was obtained. The design speed of the Rotor 67 is 16,043 rpm, the total inlet pressure is standard atmospheric pressure, the total inlet temperature is 288.15K, the design point backpressure is 122,000 Pa, the near-stall point backpressure is 126,000 Pa, and the near-blockage point backpressure is 90,000 Pa. Flutter analysis was performed for each operating condition.

[0147] Step 2: Perform structural statics and dynamics analysis to obtain the blade vibration mode shape and the natural angular frequency Specifically: is the mass matrix, is the damping matrix, is the stiffness matrix, f is the aerodynamic force, is the displacement, then the structural dynamics equation based on blade vibration is:

[0148] ,

[0149] Where: is the first derivative of displacement, which is velocity; is the second derivative of displacement, which is acceleration;

[0150] In modal analysis, aerodynamic forces are not considered. f With the damping matrix C , we get formula 1:

[0151] ,

[0152] make Formula 2:

[0153] ,

[0154] Where: is the natural angular frequency of the blade, is time, and formula 2 is substituted into formula 1 to obtain formula 3:

[0155] ,

[0156] Formula 3 is the characteristic equation. Solving Formula 3 can obtain the modal vibration shape of the blade and the natural angular frequency , and the mode shapes The corresponding modal stress and modal displacement The first-order vibration mode of the blade is obtained by calculation, and its corresponding natural frequency is f 537.9Hz;

[0157] Step 3: The calculated modal vibration shape Interpolation is applied to the CFD mesh of the blade surface, and the maximum allowable alternating stress that satisfies the high cycle fatigue life is determined for the vibration mode that is prone to flutter. , the maximum allowable vibration amplitude can be calculated Specifically:

[0158] ,

[0159] For a given design life, the maximum value of the alternating stress is taken according to the Goodman diagram of Rotor 67. for , maximum permissible vibration amplitude is the maximum value of alternating stress Divided by the modal stress Multiply by the modal displacement , modal displacement , modal stress Respectively as Figure 3 、 Figure 4 As shown, this mode is a first-order bending mode. From the figure, the maximum modal displacement is 7.04356 and the maximum modal stress is , thus obtaining the maximum allowable vibration amplitude of the blade for .

[0160] Step 4: Make the blade vibrate in a modal shape , natural angular frequency and maximum permissible vibration amplitude Perform simple harmonic oscillation and use the harmonic balance method to calculate the unsteady flow caused by blade vibration. In the unsteady flow calculation, it is necessary to consider the flow nonlinearity caused by the large amplitude of blade vibration, that is, in the unsteady flow calculation, it is necessary to analyze the frequency multiples of the blade vibration frequency, specifically:

[0161] The harmonic balance method is used to analyze the unsteady flow caused by blade vibration. The harmonic balance method can be expressed as:

[0162] ,

[0163] Where: is the variable of the unsteady flow field around the blade, is the time-averaged quantity of the unsteady flow field, For time, The harmonic number used to analyze the unsteady flow field is: i is the order corresponding to the different harmonics, and For the i The Fourier coefficients corresponding to the subharmonics are, It is i The angular frequency of the unsteady flow field corresponding to the subharmonic; the harmonic balance method is used to express the flow field variables of the unsteady flow as the time-averaged quantity plus the pulsating quantity, and approximate the pulsating quantity with the Fourier series;

[0164] The fundamental frequency of the unsteady flow caused by blade vibration is the natural frequency of the blade:

[0165] ,

[0166] In addition, the interleaf phase angle needs to be specified , interleaf phase angle Characterizes the phase difference when each blade vibrates in the traveling wave mode, that is:

[0167] ,

[0168] Where: is the vibration node diameter, is the number of leaves. When it is an even number, The value is [- / 2, / 2), where - / 2 and / 2 difference , only one of them needs to be considered in the calculation; when When is an odd number, The value is [-( -1) / 2, ( -1) / 2].

[0169] The number of blades in the entire ring of this example is 22, so its vibration node diameter The range is [-11,11). The nonlinearity caused by blade vibration is strong, and the introduction of high-order harmonics needs to be considered during calculation. =1 Harmonic analysis of the unsteady flow caused by vibration, and considering the zero node diameter, that is, the inter-blade phase angle When is 0, the input harmonic information is , 0°.

[0170] When performing unsteady calculations, if the nonlinearity of the flow field is strong, it will cause large mesh deformation and harmonic number independence verification is required.

[0171] Step 5: Based on the unsteady calculation results, calculate the work done on the blade by the unsteady aerodynamic force caused by the blade vibration during one vibration cycle, i.e., the accumulated work;

[0172] The positive or negative accumulated work value determines whether the blade meets the fatigue life requirements of the mode: if the accumulated work value is positive, it means that the work done by the airflow on the blade increases during the cycle, the amplitude will continue to increase, and there is a risk of amplitude exceeding the limit; if the accumulated work value is negative, there is no risk of amplitude exceeding the limit; specifically:

[0173] Accumulated Merit It refers to the work done on the blade by the unsteady aerodynamic force caused by blade vibration within one vibration cycle, which is:

[0174] ,

[0175] When the vibration node diameter is 0, the calculated accumulated work value is -0.2093035J, and the accumulated work density distribution on the blade surface is as follows: Figure 5 As shown in the figure, the aerodynamic accumulated work on the blade surface is not positive, so when the blade vibrates in the first-order mode and the vibration node diameter is 0, the blade amplitude will not continue to increase, meeting the fatigue life requirements.

[0176] Step 6: Traverse all possible vibration modes and inter-blade phase angles to determine whether the blade meets safety requirements.

[0177] The vibration pitch range of Rotor 67 is [-11, 11). All possible vibration modes and inter-blade phase angles are traversed, and whether the blade meets the safety requirements is determined based on the positive or negative accumulated work.

[0178] The above steps can be used to obtain the maximum allowable vibration amplitude of the blade based on the fatigue life target. This can be used for nonlinear analysis to determine whether there is a risk of blade vibration exceeding the amplitude limit, thus solving the problem of general flutter linear analysis being overly conservative.

[0179] This method, used in the numerical analysis of aerodynamic damping, subjects the blade to simple harmonic vibration at its maximum allowable amplitude. Harmonic balance methods are then used to resolve the unsteady flow caused by blade vibration. The aerodynamic accumulated work is calculated, taking into account the flow nonlinearity caused by excessive blade vibration amplitude. The positive or negative sign of the aerodynamic accumulated work is used to determine whether the blade flutter characteristics meet design requirements, completing the nonlinear analysis of blade flutter.

[0180] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention should be included in the scope of protection of the present invention.

Claims

1. A nonlinear blade flutter assessment method based on accumulated work measurement and fatigue life, characterized in that: The following steps are involved: Step 1: Establish the fluid domain calculation grid of the blade, perform steady-state calculation of the blade flow field, obtain the characteristic lines of flow rate and pressure ratio, flow rate and efficiency at different blade speeds, and thus determine the calculation conditions of blade flutter; at the same time, based on structural statics and dynamics analysis, obtain the modal vibration shapes corresponding to at least the first three vibration modes of the blade vibration and the natural angular frequency ; Step 2: Any first-order mode shape Interpolate to the calculation grid of the fluid domain of the blade, and calculate the maximum allowable vibration amplitude of the blade based on the maximum allowable alternating stress that meets the high-cycle fatigue life of the blade , specifically: The maximum allowable alternating stress that satisfies the high cycle fatigue life of the blade is obtained from the Goodman diagram , and then, the maximum allowable vibration amplitude of the blade is obtained L for: ， Where, is the modal stress, is the modal displacement; Step 3: Under the calculation conditions, let the blade vibrate in any mode , natural angular frequency and maximum permissible vibration amplitude Perform simple harmonic oscillation; Under the premise of considering the blade nonlinearity caused by the excessive amplitude in the simple harmonic vibration of the blade, the fundamental frequency and harmonic information of the unsteady flow caused by the simple harmonic vibration of the blade are specified, and the unsteady blade flutter calculation is initialized based on the steady calculation results of the blade flow field. Specifically: Specify the fundamental frequency of the unsteady flow caused by blade vibration, which is the natural frequency of the blade : ， Where, is the natural angular frequency; At the same time, specify the inter-leaf phase angle , interleaf phase angle That is, the phase difference between the blades when vibrating in the traveling wave mode, which can be expressed as: ， Where: is the vibration pitch diameter of the blade, is the number of leaves; when When it is an even number, The value is [- / 2, / 2), where - / 2 and / 2 difference , only one of them is taken during calculation; when When is an odd number, The value is [-( -1) / 2, ( -1) / 2]; In addition, to resolve the unsteady flow caused by blade vibration, the specified harmonic information is: ， in, is the harmonic number, and is a positive integer greater than 1; Therefore, the harmonic balance method is used to analyze the unsteady flow caused by blade vibration, and the unsteady pressure on the blade surface is obtained. and blade vibration speed ; Finally, the unsteady pressure on the blade surface during one vibration cycle is and blade vibration speed Integrate to get the accumulated work; Step 4: Repeat steps 2 to 3 to calculate the accumulated work corresponding to each vibration mode and evaluate the blade flutter. Specifically: If the accumulated work is positive, the amount of work done by the airflow on the blade increases during the cycle, and the amplitude will continue to increase. There is a risk of exceeding the amplitude limit, and the fatigue life requirements of the corresponding vibration mode will not be met; If the accumulated work is negative, it will not exceed the maximum allowable vibration amplitude, so there is no risk of exceeding the amplitude limit, and the fatigue life requirements of the corresponding vibration mode are met.

2. The nonlinear blade flutter assessment method based on accumulated work measurement and fatigue life according to claim 1, characterized in that: In the steady-state calculation of the blade flow field, the FANS equation for determining the calculation conditions of the blade flutter is solved, and the FANS equation is expressed in the cylindrical coordinate system as: Where: is physical time, is the conserved variable within the unit volume of the flow field, The expression is: ， is the fluid density, 、 、 are the velocity components of the fluid in the axial, circumferential and radial directions, is the total energy per unit mass of fluid; are the axial, circumferential and radial coordinate positions respectively, are the velocity components of the fluid domain calculation grid in the axial, circumferential and radial directions, 、 、 They represent the axial, circumferential, and radial fluxes generated by the movement of the computational grid in the fluid domain; are the convective fluxes of the fluid in the axial, circumferential and radial directions respectively. The convective flux represents the physical quantity transported through a unit area in unit time. The expression formula is: ， is the total enthalpy per unit mass of fluid, is the static pressure of the fluid; are the viscous fluxes of the fluid in the axial, circumferential and radial directions respectively. The viscous flux is used to describe the diffusion transport of momentum and energy caused by viscous stress and heat conduction. The expression is: ， Where, is the normal stress between fluids, is the shear stress between fluids, is the heat flow between fluids, 、 、 The heat flow between the fluids is Components in the axial, circumferential, and radial directions; is the source term, which is used to describe the influence of the internal or external flow system on the conserved quantity. The expression is: ， After solving the FANS equation, the characteristic lines of flow rate and pressure ratio, flow rate and efficiency at different blade speeds are obtained, and the calculation conditions are selected according to the characteristic lines for flutter analysis; The calculation conditions include design conditions, congestion conditions and near-stall conditions.

3. The nonlinear blade flutter assessment method based on accumulated work measurement and fatigue life according to claim 1, characterized in that: The modal vibration modes corresponding to at least the first three vibration modes of the blade vibration are obtained in step 1. and the natural angular frequency , specifically: Structural dynamics equations based on blade vibration: ， Where: is the mass matrix, is the damping matrix, is the stiffness matrix, f is the aerodynamic force acting on the blade surface, is the vibration displacement of the blade, for The first derivative of is the vibration speed of the blade; is the acceleration of the blade; In modal analysis, aerodynamic forces are not considered. f With the damping matrix C, we get formula 1: ， make Formula 2: ， Where: is the natural angular frequency of the blade, t is the physical time, and after substituting Formula 2 into Formula 1, we get Formula 3: ， Formula 3 is the characteristic equation. Solving Formula 3 will give the modal vibration shapes corresponding to each vibration mode of the blade. and the natural angular frequency , and the mode shapes The corresponding modal stress and modal displacement .

4. The nonlinear blade flutter assessment method based on accumulated work measurement and fatigue life according to claim 1, characterized in that: In step 3, the unsteady blade flutter calculation is initialized based on the steady calculation results of the blade flow field, and then the harmonic balance method is used to analyze the unsteady flow of the blade caused by the blade vibration to obtain the unsteady pressure on the blade surface corresponding to the unsteady flow of the blade. and blade vibration speed , wherein the harmonic balance method formula is: ， Where: is the variable of the unsteady flow field around the blade, is the time-averaged quantity of the unsteady flow field, is the physical time, n is the harmonic number used to analyze the unsteady flow field, i is the order corresponding to the different harmonics, and For the i The Fourier coefficients corresponding to the subharmonics are, It is i The angular frequency of the unsteady flow field corresponding to the subharmonic; the harmonic balance method is used to express the flow field variables of the unsteady flow as the time-averaged quantity plus the pulsating quantity, and approximate the pulsating quantity with the Fourier series; Harmonic balance method is used to represent the unsteady pressure on the blade surface The expression formula is: ， Where: is the time-averaged unsteady pressure on the blade surface; Similarly, the vibration speed of the blade is obtained The expression is: ， Where: is the time-averaged value of the blade vibration velocity.

5. The nonlinear blade flutter assessment method based on accumulated work measurement and fatigue life according to claim 1, characterized in that: The accumulated work described in step 3 for: ， Where, is the vibration speed of the blade, is the fluid domain mesh normal vector on the blade surface, is the differential of the leaf surface area, is the blade vibration period, dt is the physical time differential, and S is the entire blade surface area.

6. The nonlinear blade flutter assessment method based on accumulated work measurement and fatigue life according to any one of claims 1 to 5, characterized in that: The above analysis based on structural statics and dynamics obtains at least the first three modal vibration modes that are prone to flutter. and the natural angular frequency .

Images