Aerodynamic-aeroelastic coupled optimization method for bladed disk of turbomachinery considering structural factors

Abstract

The application discloses a turbomachinery blade aero-aeolastic coupling concomitant optimization method, comprising the following steps: step one, dividing the turbomachinery blade into a fluid grid to obtain the unsteady flow field information inside the turbomachinery; step two, obtaining an aero-aeolastic coupling target function; step three, obtaining design parameters; step four, establishing a concomitant harmonic balance equation to obtain concomitant variables by solving the concomitant harmonic balance equation; step five, calculating aero-aeolastic coupling sensitivity information; step six, updating the blade geometry; step seven, dividing the blade geometry into a structure grid, and carrying out finite element analysis and modal interpolation to obtain vibration modal and vibration frequency information; and step eight, when the optimized time-averaged mass flow, time-averaged total pressure ratio, time-averaged isentropic efficiency and aero-damping meet the indexes, the corresponding blade geometry is output. The application realizes the improvement of the aeolastic performance of the turbomachinery blade from both the structure and the aero-dynamics.

Description

Technical Field

[0001] This invention relates to the field of fluid-structure interaction technology for turbine blades such as aero-engines, steam turbines, and gas turbines, and specifically to an aerodynamic-aeroelastic coupling optimization method for turbine blades that incorporates structural factors. Background Technology

[0002] Flow-induced vibration in turbojet engines is a flow phenomenon resulting from the interaction between fluid and solid. Based on the difference in excitation force, it can be classified into flutter, forced vibration, and asynchronous vibration. Regardless of the type, flow-induced vibration can lead to high-cycle fatigue fracture of turbojet blades, affecting the service life of aero-engines and causing significant economic and property losses.

[0003] Due to blade vibration and static-rotation interference, the flow inside the turbomachine has inherent unsteady characteristics. To carry out research on the aerodynamic-aeroelastic coupling optimization of turbomachine blades, it is necessary to solve the control equations of the unsteady flow inside the aeroelastic system.

[0004] Currently, publicly available literature on aerodynamic-aeroelastic coupling optimization of turbomachinery blades based on the adjoint harmonic balance method can be mainly divided into two types: flutter-free optimization design for flutter instability problems, such as the paper titled "Multi-Objective Aerodynamic and Aeroelastic Coupled Design Optimization Using a Full-Viscosity Discrete Adjoint Harmonic Balance Method," and forced response minimization research for forced resonance problems, such as the paper titled "Efficient Forced Response Minimization Using a Full-Viscosity Discrete Adjoint Harmonic Balance Method." These studies can improve the aeroelastic performance of turbomachinery blades while ensuring that aerodynamic performance is not diminished or even improved.

[0005] However, to simplify the optimization process, it is assumed that the structural information of the turbine blades, such as vibration modes and frequencies, remains unchanged throughout the optimization design process. For large deformation optimization problems, this assumption will lead to a decrease in the reliability of the optimization results, requiring a re-evaluation of the results. This will increase the computational time and reduce the optimization efficiency of the aerodynamic-aeroelastic coupling optimization method for turbine blades. Summary of the Invention

[0006] To overcome the shortcomings of the existing technology, this invention provides an aerodynamic-aeroelastic coupling optimization method for turbomachine blades that integrates structural factors. This method integrates structural factors such as vibration modes and vibration frequencies into the aerodynamic-aeroelastic coupling optimization system for turbomachine blades, thereby simultaneously improving the aeroelastic performance of turbomachine blades from both structural and aerodynamic perspectives. It features high optimization efficiency, wide applicability, and high reliability of results.

[0007] To achieve the above objectives, the technical solution adopted by the present invention is as follows: An adjoint optimization method for turbine blades that integrates structural factors includes the following steps; Step 1: Divide the impeller blades into a fluid mesh, solve the unsteady harmonic equilibrium equations, and obtain the unsteady flow field information inside the impeller; Step 2: Based on the unsteady flow field information, calculate the aerodynamic performance parameters and aeroelastic performance parameters of the turbine blades respectively, and calculate the weighting coefficients of the aerodynamic-aeroelastic coupling objective function through an adaptive weighting method, and then obtain the aerodynamic-aeroelastic coupling objective function through the penalty function method; Step 3: Extract the leaf shape of the three-dimensional blade at different leaf heights and the mid-arc line and thickness distribution of the two-dimensional leaf shape at each leaf height. Keep the mid-arc line distribution unchanged, and use two sets of Hicks-Henne type functions to parameterize the perturbation of the thickness distribution in the leading and trailing edges and the middle region of the mid-arc line to obtain the design parameters. Step 4: Based on the unsteady flow field information and the aero-elastic coupling objective function, establish the adjoint harmonic balance equation based on the adjoint principle, and solve the adjoint harmonic balance equation to obtain the adjoint variables; Step 5: Calculate the aerodynamic-aeroelastic coupling sensitivity information based on the aerodynamic-aeroelastic coupling objective function, the design parameters, and the accompanying variables; Step 6: Update the design parameters based on the sequential quadratic programming method and the aerodynamic-aeroelastic coupling sensitivity information, and update the blade geometry using the updated design parameter information through the Hicks-Henne type function method; Step 7: When the turbine blades meet the structural information update criteria, perform structural mesh generation on the blade geometry, and carry out finite element analysis and modal interpolation to obtain its vibration modes and vibration frequency information; Step 8: Repeat steps 1 to 7. When the optimized time-averaged mass flow rate, time-averaged total pressure ratio, and time-averaged isentropic efficiency meet the optimization criteria, output the corresponding blade geometry.

[0008] Furthermore, the fluid mesh generation in step one is performed using the AUTOGRID5 module of the software NUMECA; The specific expression for the unsteady harmonic equilibrium equation is as follows:

[0009] In the formula, Q represents the flow field variable; R represents the spatial residual; and E represents the time spectrum source term. The obtained flow field variables Q include pressure p, temperature T, and density. And velocity v, which is the unsteady flow field information.

[0010] Furthermore, the aerodynamic performance parameters in step two include mass flow rate. Pressure ratio and efficiency The aeroelastic performance parameters are the accumulated work. The specific calculation formula is as follows:

[0011] In the formula, ds is the area of ​​the control unit; For unsteady aerodynamic loads on the blade surface; For the grid movement speed; The unit outward direction of the blade surface; The specific implementation steps of the adaptive weighting method are as follows: First, the dimensionless time-averaged values ​​of the blade aerodynamic performance parameters are calculated, and the specific expressions are as follows:

[0012]

[0013]

[0014] In the formula, The time-average entropy efficiency is dimensionless; The time-averaged mass flow rate is dimensionless; The dimensionless time-averaged total voltage ratio; N is the number of harmonics; This is the initial value for isentropic efficiency; This is the initial value of the total pressure ratio; This is the initial value for the mass flow rate; Let i be the i-th time point; Secondly, the aerodynamic damping, used to characterize the aeroelastic performance parameters, is calculated, with the specific expression as follows:

[0015] In the formula, Modal displacement; ω is the angular frequency of the blade vibration; To accumulate merit; For aerodynamic damping; This is the initial aerodynamic damping; Pi; Finally, the weighting coefficients for aerodynamic and aeroelastic properties are determined based on the dimensionless time-averaged isentropic efficiency and aerodynamic damping, as shown in the following expression:

[0016]

[0017] In the formula, The weighting coefficient for the dimensionless time-averaged isentropic efficiency; The weighting factor for aerodynamic damping; The penalty function method uses a quadratic function to constrain the time-averaged mass flow rate and the time-averaged total pressure ratio, and the specific expression is as follows:

[0018]

[0019] In the formula, For time-averaged mass flow rate constraints; The time-averaged total pressure ratio is a constraint. The time-average mass flow rate penalty function coefficient; The time-averaged total pressure ratio penalty function coefficient; The specific expression of the aerodynamic-aeroelastic coupling objective function is as follows:

[0020] In the formula, The objective function is the aerodynamic-aeroelastic coupling function.

[0021] Furthermore, in step three, the leaf shape at different leaf heights and the mid-arc line and thickness distribution of each two-dimensional leaf shape at each leaf height are extracted, as specifically expressed below:

[0022] In the formula, the superscript k represents different leaf heights; the superscript chord represents the middle arc. Represents the circumferential coordinates of the two-dimensional leaf shape at the k-th leaf height; Represents the mid-arc coordinates of the two-dimensional leaf profile at the k-th leaf height; This represents the thickness distribution of the two-dimensional leaf shape at the k-th leaf height; Keep the distribution of the middle arc unchanged, that is, keep Without changing the equation, the thickness distribution of the two-dimensional airfoil can be further expressed as follows:

[0023] In the formula, This represents the blade thickness distribution before optimization. The disturbance amount is the thickness distribution of the leading and trailing edges of the two-dimensional airfoil; This refers to the disturbance amount in the central region of the two-dimensional airfoil. In step three, the leading and trailing edges of the blade are parameterized based on the Hicks-Henne type function method, and the specific expressions are as follows:

[0024] In the formula, This is the j-th design parameter for the k-th leaf height; The perturbation amount is the thickness distribution of the leading or trailing edge of the two-dimensional leaf shape at the k-th leaf height; Let be the dimensionless radial coordinate of the j-th design parameter; Let be the dimensionless axial coordinate of the j-th design parameter; For dimensionless radial coordinates; For dimensionless axial coordinates; subscript j indicates design parameter number; c indicates coordinate center; le indicates leading edge; te indicates trailing edge; d indicates dimensionless. , , and The definition is as follows:

[0025]

[0026]

[0027]

[0028] In the formula, the subscripts h and s represent the blade root and blade tip of the turbomachinery blade, respectively; M and s represent the number of different blade heights and the number of design parameters for the same blade height, respectively. The parameterization of the central region of the blade's arc using the Hicks-Henne type function method is as follows:

[0029] In the formula, This is the j-th design parameter for the k-th leaf height; The disturbance amount of the thickness distribution of the j-th design parameter in the middle region of the two-dimensional airfoil with the k-th airfoil height; and The definition is as follows:

[0030]

[0031] From the two sets of Hicks-Henne type functions mentioned above, we can see that the total number of design parameters is s+2.

[0032] The specific steps for establishing the accompanying harmonic balance equation in step four are as follows: First, the aerodynamic-aeroelastic coupling optimization of the turbomachinery blades is written in the following notation form: Objective function:

[0033] constraint:

[0034] Secondly, the constraints are transformed into the objective function using the Lagrange multiplier method, as shown in the following expression:

[0035] In the formula, J is the Lagrange function. These are Lagrange multipliers, also known as adjoint variables; Finally, the partial derivative of J with respect to the flow field variable Q is calculated, and the specific expression is as follows:

[0036] By rearranging and combining like terms, we can obtain the expression for the unsteady adjoint equation as follows:

[0037] By iteratively solving the above adjoint equation, the adjoint variables can be obtained. .

[0038] The calculation of aerodynamic-aeroelastic coupling sensitivity information in step five is achieved by calculating J with respect to... The derivative is obtained, and the specific expression is as follows: .

[0039] In the formula, This is the total design parameter vector; This refers to the aerodynamic-aeroelastic coupling sensitivity information.

[0040] The sequential quadratic programming method used in step six to update the design parameters is specifically calculated using the following formula:

[0041] In the formula, is the partial derivative of the objective function I with respect to the design parameters; h is the constraint of the flow field equation; The flow field equations with respect to the design parameters The partial derivatives; is the multiplier of the Lagrange function; the superscript T indicates transpose; the subscript k indicates the iteration step; h The specific expression is as follows:

[0042]

[0043]

[0044] The specific implementation process for updating the blade geometry is as follows: First, the perturbation amount of the mid-curve distribution is calculated from the updated design parameters; then, the perturbation amount is superimposed on the mid-curve distribution of the original airfoil to obtain the mid-curve distribution of the optimized airfoil; finally, the thickness distribution of the original airfoil is superimposed on the mid-curve distribution of the optimized airfoil to update the blade geometry.

[0045] The structural update criterion in step seven is based on the L2 norm of the design parameters, and the specific expression is as follows:

[0046] In the formula, The design parameter with the largest value is considered to have a large blade profile change when its L2 norm is greater than the threshold of the largest design parameter (0.1 times). Otherwise, the blade profile change is considered small and no structural information needs to be updated.

[0047] The structural mesh was generated using the Workbench module of ANSYS software, with the sweep method selected and hexahedral mesh elements chosen. The modal interpolation is performed using the radial basis function method to interpolate the vibration modes from the solid mesh in step seven to the fluid mesh in step one. The specific expression is as follows:

[0048] In the formula, m is the number of interpolation nodes; is the interpolation node, and is the fluid mesh coordinate in step one; x is the displacement vector of the original node, and is the solid mesh coordinate in step seven; The weighting coefficients for interpolation node i; These are the interpolation basis functions; These are the interpolated vibration modes.

[0049] The specific expressions for the time-averaged mass flow rate, time-averaged total pressure ratio, time-averaged isentropic efficiency, and aerodynamic damping in step eight are as follows:

[0050]

[0051]

[0052]

[0053] In the formula, the subscript opt ​​represents the optimized performance parameters; org represents the unoptimized performance parameters. The optimized time-averaged mass flow rate; The average mass flow rate before optimization; The optimized average total pressure ratio; The average total pressure ratio before optimization; This refers to the optimized time-averaged isentropic efficiency; The time-average isentropic efficiency before optimization; It is for aerodynamic damping.

[0054] The beneficial effects of this invention are: This invention considers the influence of structural factors such as vibration modes and vibration frequencies on the basis of an aerodynamic-aeroelastic coupled optimization system. On the one hand, it can ensure the reliability of optimization results and improve the applicability of the optimization system for large deformation optimization problems; on the other hand, it can improve the aeroelastic performance of turbine blades from both aerodynamic and structural perspectives, thereby significantly enhancing the optimization effect.

[0055] This invention proposes an update criterion for blade structure information. This criterion determines the magnitude of blade shape change based on the ratio of the L2 norm of the design parameters to the maximum design parameter. The value of this ratio relative to 0.1 determines whether to call ANSYS for finite element analysis. The advantage is that it reduces the additional computation time caused by calling ANSYS for modal analysis in each optimization step, making it more efficient.

[0056] This invention proposes an adaptive weighted aerodynamic-aeroelastic sensitivity coupling strategy for turbojet blades. This strategy automatically adjusts the weighting coefficients based on the ratio of aerodynamic and aeroelastic performance parameters according to the optimization objective, simplifying the calculation of the aerodynamic-aeroelastic coupling objective function.

[0057] This invention employs a Python script to integrate the following steps: solving the unsteady harmonic balance equation in step one, solving the aerodynamic-aeroelastic coupling objective function in step two, parameterizing the turbine blades based on the Hickss-Henne type function in step three, solving the adjoint harmonic balance equation in step four, analyzing the aerodynamic-aeroelastic coupling sensitivity information based on the adjoint method in step five, updating the design parameters based on the sequential quadratic programming method in step six, and updating the structural information in step seven. This achieves full automation of aerodynamic-aeroelastic coupling optimization, improves optimization efficiency, and reduces the reliance on the designer's experience in the optimization process. Attached Figure Description

[0058] Figure 1 This is a schematic diagram of the optimized process of the present invention; Figure 2 A schematic diagram illustrating the specific process of blade parameterization; Figure 3A schematic diagram illustrating the convergence history of the optimization of aerodynamic performance-isentropic efficiency and aeroelastic performance-accumulated work; Figure 4 This is a schematic diagram comparing the original leaf shape and the optimized leaf shape at 50% leaf height. Detailed Implementation

[0059] The present invention will now be described in further detail with reference to the accompanying drawings.

[0060] like Figure 1 As shown, this invention discloses a method for aerodynamic-aeroelastic coupling optimization of turbofan blades that incorporates structural factors, comprising the following steps: Step 1: Use the AUTOGRID5 module of NUMECA software to generate a fluid mesh for the turbine blades and solve the unsteady harmonic equilibrium equations to obtain the unsteady flow field information inside the turbine. The expression for the solved unsteady harmonic equilibrium equations is as follows:

[0061] In the formula, Q represents the flow field variable; R represents the spatial residual; and E represents the time spectrum source term. The obtained flow field variables Q include pressure p, temperature T, and density. Information on unsteady flow fields with velocity v.

[0062] Step 2: Based on the unsteady flow field information obtained in Step 1, calculate the mass flow rate of the turbine blades. Pressure ratio and efficiency Aerodynamic performance parameters and accumulated work, etc., are considered. The weighting coefficients of the aerodynamic-aeroelastic coupling objective function are calculated using an adaptive weighting method. Finally, the aerodynamic-aeroelastic coupling objective function is obtained using the penalty function method. The formula for calculating the aeroelastic performance parameter - accumulated work is:

[0063] In the formula, ds is the area of ​​the control unit; For unsteady aerodynamic loads on the blade surface; For the grid movement speed; The unit outward direction of the blade surface; The specific implementation steps of the adaptive weighting method are as follows: First, the dimensionless time-averaged values ​​of the blade aerodynamic performance parameters are calculated, and the specific expressions are as follows:

[0064]

[0065]

[0066] In the formula, The time-average entropy efficiency is dimensionless; The time-averaged mass flow rate is dimensionless; The dimensionless time-averaged total voltage ratio; N is the number of harmonics; This is the initial value for isentropic efficiency; This is the initial value of the total pressure ratio; This is the initial value for the mass flow rate; Let i be the i-th time point; Secondly, the aerodynamic damping used to characterize the aeroelastic performance is calculated, with the specific expression as follows:

[0067] In the formula, Modal displacement; ω is the angular frequency of the blade vibration; To accumulate merit; For aerodynamic damping; This is the initial aerodynamic damping; Pi; Finally, the weighting coefficients for aerodynamic and aeroelastic properties are determined based on the dimensionless time-averaged efficiency and aerodynamic damping, as shown in the following expression:

[0068]

[0069] In the formula, The weighting coefficients for the dimensionless time-averaged efficiency; The weighting factor for aerodynamic damping; The penalty function method uses a quadratic function to constrain the mass flow rate and pressure ratio, and the specific expression is as follows:

[0070]

[0071] In the formula, For mass flow constraints; The time-averaged total pressure ratio is a constraint. The time-average mass flow rate penalty function coefficient; The time-averaged total pressure ratio penalty function coefficient; The specific expression for the aerodynamic-aeroelastic coupling objective function is shown below:

[0072] In the formula, The objective function is the aerodynamic-aeroelastic coupling function. Step 3: Extract the leaf shape of the three-dimensional blade at different leaf heights and the mid-curve and thickness distribution of the two-dimensional leaf shape at each leaf height. Keeping the thickness distribution unchanged, two sets of Hicks-Henne type functions are used to parameterize the perturbation amounts of the leading and trailing edges and the middle region of the mid-curve, respectively, to obtain the design parameters. The trailing and leading edges of the blade are parameterized based on the Hicks-Henne type function method, with the specific expressions as follows:

[0073] In the formula, This is the j-th design parameter for the k-th leaf height; denoted as circumferential disturbance at the leading or trailing edge of the k-th blade height; r is the radial coordinate; x is the axial coordinate; subscript j indicates the design parameter number; c indicates the coordinate center; le indicates the leading edge; te indicates the trailing edge; d indicates dimensionless design. The parameterization of the middle part of the blade's arc is based on the Hicks-Henne type function method, and the specific expression is as follows:

[0074] In the formula, This is the j-th design parameter for the k-th leaf height; Let be the disturbance amount of the thickness distribution of the j-th design parameter in the central region of the two-dimensional airfoil with the k-th blade height; s is the number of design parameters; m and n are defined as follows:

[0075]

[0076] From the two sets of Hicks-Henne type functions mentioned above, we can see that the total number of design parameters is s+2; Step 4: Based on the unsteady flow field information obtained in Step 1, the aerodynamic-aeroelastic coupling objective function obtained in Step 2, and the adjoint principle, establish the adjoint harmonic balance equation and solve for the adjoint variables. The specific steps for establishing the adjoint harmonic balance equation are as follows: First, the aerodynamic-aeroelastic coupling optimization of the turbomachinery blades is written in the following notation form: Objective function:

[0077] constraint:

[0078] Secondly, the constraints are transformed into the objective function using the Lagrange multiplier method, as shown in the following expression:

[0079] In the formula, These are Lagrange multipliers, also known as adjoint variables; Finally, the partial derivative of J with respect to the flow field variable Q is calculated, and the specific expression is as follows:

[0080] By rearranging and combining like terms, we can obtain the expression for the unsteady adjoint equation as follows:

[0081] By iteratively solving the above adjoint equation, the adjoint variables can be obtained. ; Step 5: Based on the aerodynamic-aeroelastic coupling objective function obtained in Step 2, the design parameters obtained in Step 3, and the adjoint variables obtained in Step 4, calculate the aerodynamic-aeroelastic coupling sensitivity information. The specific calculation formula for the sequential quadratic programming method is as follows:

[0082] In the formula, is the partial derivative of the objective function I with respect to the design parameters; h is the constraint of the flow field equation; The flow field equations with respect to the design parameters The partial derivatives; is the multiplier of the Lagrange function; the superscript T indicates transpose; the subscript k indicates the iteration step; h The specific expression is as follows:

[0083]

[0084]

[0085] The specific implementation process for updating the blade geometry is as follows: First, the disturbance amount of the mid-curve distribution is calculated based on the updated design parameters; then, the disturbance amount is superimposed on the mid-curve distribution of the original blade profile to obtain the mid-curve distribution of the optimized blade profile; finally, the thickness distribution of the original blade profile is superimposed on the mid-curve distribution of the optimized blade profile to update the blade geometry. Step 7: Once the turbine blades meet the structural information update criteria, the blade geometry is meshed using the Workbench module of ANSYS software. The meshing method is sweep, and hexahedral mesh elements are selected. Finite element analysis and modal interpolation are then performed to obtain the vibration modes and frequencies. The structural update criteria are based on the L2 norm of the design parameters, and the specific expression is shown below:

[0086] In the formula, The design parameter with the largest value; Modal interpolation is implemented using the radial basis function method, with the specific expression as follows:

[0087] In the formula, m is the number of interpolation nodes; is the interpolation node, and is the fluid mesh coordinate in step one; x is the displacement vector of the original node, and is the solid mesh coordinate in step seven; The weighting coefficients for interpolation node i; These are the interpolation basis functions; These are the interpolated vibration modes; Step 8: Repeat steps 1 to 7. When the optimized time-averaged mass flow rate, time-averaged total pressure ratio, time-averaged isentropic efficiency, and aerodynamic damping meet the optimization targets, output the corresponding blade geometry. The specific expression is as follows:

[0088]

[0089]

[0090]

[0091] In the formula, the subscript opt ​​represents the optimized performance parameters; org represents the unoptimized performance parameters.

[0092] To illustrate the purpose, technical solution, and advantages of this invention in more detail, the NASA Rotor 67 transonic fan rotor is used as an example for further detailed explanation: Example: like Figure 2 , Figure 3 , Figure 4 As shown, the flow field solver and adjoint solver in this example are implemented based on self-written Fortran programs TurboXD and TurboADJ, and the optimization platform is implemented based on a self-written Python script. Finite element mesh generation and analysis are implemented using ANSYS Workbench, the modal interpolation program is implemented based on a self-written Python script, blade parameterization is implemented using a self-written Fortran program, and the mesh generation program is NUMECA's AutoGrid5. The specific implementation steps of this example are as follows: S01. The fluid computational domain was meshed using NUMECA's AutoGrid5, and unsteady aeroelastic flow field analysis was performed based on the self-developed flow field solver - TurboXD. The unsteady flow field information was obtained and the flow field data was saved in solution.save. Based on the unsteady flow field information, the aeroelastic performance parameters - accumulated work and aerodynamic performance parameters - mass flow rate, pressure ratio and isentropic efficiency were calculated respectively. The aerodynamic-aeroelastic coupling objective function was calculated based on the adaptive weighting criterion. The mass flow rate and pressure ratio were constrained respectively by the penalty function method. The penalty function coefficients of mass flow rate and pressure ratio were 200. S02. Using a self-developed blade parameterization program, two-dimensional blade profiles were extracted at 10%, 20%, 35%, 50%, 65%, 80%, and 90% of the blade height in the three-dimensional blade design. Further, the mid-curve distribution and thickness distribution of the two-dimensional blade profiles were extracted. The mid-curve was divided into 10 equal parts along the axial direction (including the leading and trailing edges), and two different Hicks-Henne type function methods were used to parameterize the perturbation of the mid-curve. Figure 2 The process of blade parameterization is demonstrated. The blade parameterization information is written to the design_var_inf.dat file, and the initial values ​​of the design parameters are written to the design_var.dat file. Usually, the initial values ​​of the design parameters are set to 0. S03. Input the unsteady aeroelastic flow field information and aero-aeroelastic coupling objective function obtained in S01 into the self-compiled adjoint solver TurboADJ, and perform adjoint field analysis to obtain adjoint variables. Substitute the adjoint variables into the aero-aeroelastic coupling sensitivity calculation formula to obtain the sensitivity information of the aero-aeroelastic coupling objective function with respect to the design parameters obtained in S02, and write it into the sen_inf.dat file. S04. Using a self-written Python optimization script, read the gradient information from the sen_inf.dat file obtained in S03, update the design parameters based on the sequential quadratic programming method, write the updated design parameters into the design_var.dat file, update the perturbation amount of the arc line in the blade according to the updated design parameters, superimpose the updated design parameters onto the original arc line distribution to obtain the updated arc line distribution, and superimpose the original thickness distribution to obtain the updated blade shape. S05. Based on the optimization objective, determine whether the aerodynamic-aeroelastic coupling performance of the updated blade geometry obtained in S04 meets the performance indicators. If the performance indicators are met, output the final blade profile; otherwise, proceed with the next optimization step based on the magnitude of the blade profile change. S06. Calculate the L2 norm of the updated design parameters obtained in S04 and compare it with the design parameter with the largest absolute value. When the ratio of the L2 norm of the design parameter to the design parameter with the largest absolute value is less than 0.1, the design parameter information does not need to be updated, and S01 to S05 are repeated; when the ratio of the L2 norm of the design parameter to the design parameter with the largest absolute value is greater than 0.1, the vibration mode and vibration frequency information need to be updated based on ANSYS, and the vibration mode is interpolated from the structural mesh to the fluid mesh using the radial basis function method, and S01 to S05 are repeated. During the optimization process, the convergence history of aerodynamic performance parameters and aeroelastic performance parameters is as follows: Figure 3 As shown in the figure, after 19 optimization steps, the accumulated work decreased by about 10%, while aerodynamic performance parameters such as isentropic efficiency, mass flow rate, and pressure ratio remained basically unchanged. A comparison of the blade geometry at 50% blade height before and after optimization is also provided. Figure 4 As shown in the diagram. The solid black line represents the pre-optimized leaf shape, and the dashed line represents the post-optimized leaf shape. From... Figure 4 It can be seen that, compared to the unoptimized version, the optimized blade profile has increased camber at the leading edge and decreased camber in the middle and trailing edges. The change in blade shape indicates that shifting the blade load rearward is beneficial for improving the aerodynamic and aeroelastic performance of the fan blades.

[0093] The above description is merely 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 within the protection scope of the present invention.

Claims

1. A method for aerodynamic-aeroelastic coupling optimization of turbomachinery blades incorporating structural factors, characterized in that, Includes the following steps; Step 1: Divide the impeller blades into a fluid mesh, solve the unsteady harmonic equilibrium equations, and obtain the unsteady flow field information inside the impeller; Step 2: Based on the unsteady flow field information, calculate the aerodynamic performance parameters and aeroelastic performance parameters of the turbine blades respectively, and calculate the weighting coefficients of the aerodynamic-aeroelastic coupling objective function through an adaptive weighting method, and then obtain the aerodynamic-aeroelastic coupling objective function through the penalty function method; Step 3: Extract the leaf shape of the three-dimensional blade at different leaf heights and the mid-arc line and thickness distribution of the two-dimensional leaf shape at each leaf height. Keep the mid-arc line distribution unchanged, and use two sets of Hicks-Henne type functions to parameterize the perturbation of the thickness distribution in the leading and trailing edges and the middle region of the mid-arc line to obtain the design parameters. Step 4: Based on the unsteady flow field information and the aero-elastic coupling objective function, establish the adjoint harmonic balance equation based on the adjoint principle, and solve the adjoint harmonic balance equation to obtain the adjoint variables; Step 5: Calculate the aerodynamic-aeroelastic coupling sensitivity information based on the aerodynamic-aeroelastic coupling objective function, the design parameters, and the accompanying variables; Step 6: Update the design parameters based on the sequential quadratic programming method and the aerodynamic-aeroelastic coupling sensitivity information, and update the blade geometry using the updated design parameter information through the Hicks-Henne type function method; Step 7: When the turbine blades meet the structural information update criteria, perform structural mesh generation on the blade geometry, and carry out finite element analysis and modal interpolation to obtain its vibration modes and vibration frequency information; Step 8: Repeat steps 1 to 7. When the optimized time-averaged mass flow rate, time-averaged total pressure ratio, time-averaged isentropic efficiency, and aerodynamic damping meet the optimization criteria, output the corresponding blade geometry.

2. The aerodynamic-aeroelastic coupling optimization method for turbojet blades incorporating structural factors as described in claim 1, characterized in that, The fluid mesh generation in step one is performed using the AUTOGRID5 module of the software NUMECA. The specific expression for the unsteady harmonic equilibrium equation is as follows: In the formula, Q represents the flow field variable; R represents the spatial residual; and E represents the time spectrum source term. The obtained flow field variables Q include pressure p, temperature T, and density. And velocity v, which is the unsteady flow field information.

3. The aerodynamic-aeroelastic coupling optimization method for turbojet blades based on integrated structural factors as described in claim 2, characterized in that, The aerodynamic performance parameters in step two include mass flow rate. Pressure ratio and efficiency The aeroelastic performance parameters are the accumulated work. The specific calculation formula is as follows: In the formula, ds is the area of ​​the control unit; For unsteady aerodynamic loads on the blade surface; For the grid movement speed; The unit outward direction of the blade surface; The specific implementation steps of the adaptive weighting method are as follows: First, the dimensionless time-averaged values ​​of the blade aerodynamic performance parameters are calculated, and the specific expressions are as follows: In the formula, The time-average entropy efficiency is dimensionless; The time-averaged mass flow rate is dimensionless; The dimensionless time-averaged total voltage ratio; N is the number of harmonics; This is the initial value for isentropic efficiency; This is the initial value of the total pressure ratio; This is the initial value for the mass flow rate; Let i be the i-th time point; Secondly, the aerodynamic damping, used to characterize the aeroelastic performance parameters, is calculated, with the specific expression as follows: In the formula, Modal displacement; ω is the angular frequency of the blade vibration; To accumulate merit; For aerodynamic damping; This is the initial aerodynamic damping; Pi; Finally, the weighting coefficients for aerodynamic and aeroelastic properties are determined based on the dimensionless time-averaged isentropic efficiency and aerodynamic damping, as shown in the following expression: In the formula, The weighting coefficient for the dimensionless time-averaged isentropic efficiency; The weighting factor for aerodynamic damping; The penalty function method uses a quadratic function to constrain the time-averaged mass flow rate and the time-averaged total pressure ratio, and the specific expression is as follows: In the formula, For time-averaged mass flow rate constraints; The time-averaged total pressure ratio is a constraint. The time-average mass flow rate penalty function coefficient; The time-averaged total pressure ratio penalty function coefficient; The specific expression of the aerodynamic-aeroelastic coupling objective function is as follows: In the formula, The objective function is the aerodynamic-aeroelastic coupling function.

4. The aerodynamic-aeroelastic coupling optimization method for turbojet blades based on the fusion of structural factors as described in claim 3, characterized in that, In step three, the leaf shape at different leaf heights and the mid-arc line and thickness distribution of each two-dimensional leaf shape at each leaf height are extracted. The specific expressions are as follows: In the formula, the superscript k represents different leaf heights; The superscript chord indicates a mid-curve; Represents the circumferential coordinates of the two-dimensional leaf shape at the k-th leaf height; Represents the mid-arc coordinates of the two-dimensional leaf profile at the k-th leaf height; This represents the thickness distribution of the two-dimensional leaf shape at the k-th leaf height; Keep the distribution of the middle arc unchanged, that is, keep Without changing the equation, the thickness distribution of the two-dimensional airfoil can be further expressed as follows: In the formula, This represents the blade thickness distribution before optimization. This represents the perturbation amount of the thickness distribution at the leading and trailing edges of a two-dimensional airfoil; This represents the disturbance in the central region of a two-dimensional leaf shape. In step three, the leading and trailing edges of the blade are parameterized using the Hicks-Henne type function method, and the specific expressions are as follows: In the formula, This is the j-th design parameter for the k-th leaf height; This represents the perturbation amount of the thickness distribution of the leading or trailing edge of the two-dimensional leaf shape at the k-th leaf height; Let be the dimensionless radial coordinate of the j-th design parameter; Let be the dimensionless axial coordinate of the j-th design parameter; For dimensionless radial coordinates; For dimensionless axial coordinates; subscript j indicates design parameter number; c indicates coordinate center; le indicates leading edge; te indicates trailing edge; d indicates dimensionless. , , and The definition is as follows: In the formula, the subscripts h and s represent the blade root and blade tip of the turbomachinery blade, respectively; M and s represent the number of different blade heights and the number of design parameters for the same blade height, respectively. The parameterization of the central region of the blade's arc using the Hicks-Henne type function method is as follows: In the formula, This is the j-th design parameter for the k-th leaf height; The disturbance amount of the thickness distribution of the j-th design parameter in the middle region of the two-dimensional airfoil with the k-th airfoil height; and The definition is as follows: From the two sets of Hicks-Henne type functions mentioned above, we can see that the total number of design parameters is s+2.

5. The aerodynamic-aeroelastic coupling optimization method for turbojet blades based on integrated structural factors according to claim 4, characterized in that, The specific steps for establishing the accompanying harmonic balance equation in step four are as follows: First, the aerodynamic-aeroelastic coupling optimization of the turbomachinery blades is written in the following notation form: Objective function: constraint: Secondly, the constraints are transformed into the objective function using the Lagrange multiplier method, as shown in the following expression: In the formula, J is the Lagrange function. These are Lagrange multipliers, also known as adjoint variables; Finally, the partial derivative of J with respect to the flow field variable Q is calculated, and the specific expression is as follows: By rearranging and combining like terms, we can obtain the expression for the unsteady adjoint equation as follows: By iteratively solving the above adjoint equation, the adjoint variables can be obtained. .

6. The aerodynamic-aeroelastic coupling optimization method for turbojet blades based on integrated structural factors according to claim 5, characterized in that, The aerodynamic-aeroelastic coupling sensitivity information calculation in step five is achieved by calculating J with respect to... The derivative is obtained, and the specific expression is as follows: ； In the formula, This is the total design parameter vector; This refers to the aerodynamic-aeroelastic coupling sensitivity information.

7. The aerodynamic-aeroelastic coupling optimization method for turbojet blades incorporating structural factors according to claim 6, characterized in that, The sequential quadratic programming method used in step six to update the design parameters is specifically calculated using the following formula: In the formula, is the partial derivative of the objective function I with respect to the design parameters; h is the constraint of the flow field equation; The flow field equations with respect to the design parameters The partial derivatives; is the multiplier of the Lagrange function; the superscript T indicates transpose; the subscript k indicates the iteration step; h The specific expression is as follows: The specific implementation process for updating the blade geometry is as follows: First, the perturbation amount of the mid-curve distribution is calculated from the updated design parameters; then, the perturbation amount is superimposed on the mid-curve distribution of the original airfoil to obtain the mid-curve distribution of the optimized airfoil; finally, the thickness distribution of the original airfoil is superimposed on the mid-curve distribution of the optimized airfoil to update the blade geometry.

8. The aerodynamic-aeroelastic coupling optimization method for turbojet blades based on integrated structural factors according to claim 7, characterized in that, The structural update criterion in step seven is based on the L2 norm of the design parameters, and the specific expression is as follows: In the formula, The design parameter with the largest value is considered to have a large blade shape change when the L2 norm of the design parameter is greater than the threshold of the largest design parameter, i.e., 0.1 times. Otherwise, the blade shape change is considered to be small and the structural information does not need to be updated.

9. The aerodynamic-aeroelastic coupling optimization method for turbojet blades based on integrated structural factors as described in claim 8, characterized in that, The structural mesh was generated using the Workbench module of ANSYS software, with the sweep method selected and hexahedral mesh elements chosen. The modal interpolation uses the radial basis function method to interpolate the vibration modes from the solid mesh in step seven to the fluid mesh in step one. The specific expression is as follows: In the formula, m is the number of interpolation nodes; is the interpolation node, and is the fluid mesh coordinate in step one; x is the displacement vector of the original node, and is the solid mesh coordinate in step seven; The weighting coefficients for interpolation node i; These are the interpolation basis functions; These are the interpolated vibration modes.

10. The aerodynamic-aeroelastic coupling optimization method for turbojet blades based on fused structural factors according to claim 9, characterized in that, The specific expressions for the optimization indices of time-averaged mass flow rate, time-averaged total pressure ratio, time-averaged isentropic efficiency, and aerodynamic damping in step eight are as follows: In the formula, the subscript opt ​​represents the optimized performance parameters; org represents the unoptimized performance parameters. The optimized time-averaged mass flow rate; The average mass flow rate before optimization; The optimized average total pressure ratio; The average total pressure ratio before optimization; This refers to the optimized time-averaged isentropic efficiency; The time-average isentropic efficiency before optimization; It is for aerodynamic damping.

Images