Large-bypass-ratio fan blade multidisciplinary coupling adjoint optimization method

Abstract

The invention discloses a multidisciplinary coupling adjoint optimization method for a large-bypass-ratio fan blade. The method comprises the following steps of: 1, extracting two-dimensional blade profiles of different blade heights of a large-bypass-ratio fan blade and blade modeling parameters of each blade height; 2, obtaining optimization design parameters; 3, obtaining an unsteady harmonic balance equation according to the unsteady flow field equation, and calculating a pneumatic-aeroelastic multidisciplinary coupling objective function; 4, establishing an unsteady adjoint equation, and performing iterative solution to obtain an unsteady adjoint variable; 5, obtaining sensitivity information of the aerodynamic-aeroelastic multidisciplinary coupling objective function about design parameters; step 6, updating and optimizing design parameters; 7, updating blade modeling parameters and blade geometry, and updating the computational grid by adopting a linear elasticity method; and 8, when the two norms of the sensitivity information meet the convergence standard, outputting the optimized blade modeling parameters and blade geometry. The method has the advantages of being low in calculation cost, high in optimization design automation degree and high in optimization efficiency.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of aeroengine blade aerodynamic and aeroelastic multidisciplinary coupling, and particularly relates to a large-bypass-ratio fan blade multidisciplinary coupling concomitant optimization method. BACKGROUND

[0002] The design requirements of high load, high efficiency and light weight of aeroengines lead to the frequent occurrence of flow-induced vibration problems such as flutter and forced response of large-bypass-ratio fan blades. Flutter and forced response can cause high-cycle fatigue fracture of fan blades, and can even cause engine damage and personnel injury accidents. Therefore, it is of great importance to develop an efficient aerodynamic-aeroelastic multidisciplinary coupling optimization method for large-bypass-ratio fan blades and carry out corresponding optimization design research, so as to reduce the risk of flow-induced vibration of fan blades and improve the operation safety of aeroengines.

[0003] According to whether gradient information is needed, the large-bypass-ratio fan blade multidisciplinary coupling optimization method is mainly divided into two categories: global optimization method and local optimization method. The global optimization method has a "dimension curse" problem, that is, the calculation amount involved in optimization increases in geometric progression with the increase of the number of design variables. For multi-design-variable optimization problems, the global optimization method is relatively inefficient. Gradient optimization methods based on concomitant models are relatively efficient for multi-design-variable optimization problems because the number of flow field calculations required for sensitivity calculation is independent of the number of design variables, and therefore are favored by researchers and have been widely studied.

[0004] At present, the local optimization method based on concomitant sensitivity analysis either focuses on a single structural performance, such as the fan blade maximum equivalent stress optimization method based on the concomitant model of the patent with the publication number CN118313218A; or focuses on a single aerodynamic performance, such as the patent with the publication number CN115358167A, which is an integrated aerodynamic concomitant optimization design method for engines; or although both aerodynamic and aeroelastic coupling performances are considered, it is only for single-bypass-ratio fans, such as the papers entitled "Multi-Objective Aerodynamic and Aeroelastic Coupled Design Optimization Using a Full Viscosity Discrete Adjoint Harmonic Balance Method" and "Concurrent Blade Aerodynamic-Aero-elastic Design Optimization Using Adjoint Method".

[0005] The above research status has seriously hindered the development and progress of modern advanced multidisciplinary coupled optimization design technology for high-bypass ratio fan blades. Therefore, developing a multidisciplinary coupled optimization method suitable for high-bypass ratio fan blades and conducting corresponding optimization design research have important engineering practical value. Summary of the Invention

[0006] In order to overcome the defects of the above-mentioned existing technologies, the present invention provides a multidisciplinary coupling accompanying optimization method for large bypass ratio fan blades. This method can improve the aeroelastic stability of large bypass ratio fan blades within the full working range while ensuring that the aerodynamic performance of the internal and external bypass is not reduced. It has the characteristics of low computational cost, high degree of automation in optimization design and high optimization efficiency.

[0007] In order to achieve the above object, the technical solution adopted by the present invention is: A multidisciplinary coupled adjoint optimization method for a high bypass ratio fan blade comprises the following steps: Step 1: To maximize the design freedom and ensure that the blade shapes at the root, tip, and mid-blade are all considered, extract the two-dimensional blade shapes of at least seven different blade heights for the high bypass ratio fan blade and extract the blade shape parameters for each blade height; Step 2: Use the Hicks-Henne function method to parameterize the radial perturbation of the blade modeling parameters obtained in step 1 to obtain the optimized design parameters; Step 3: Based on the unsteady flow field equation, the unsteady flow field variables are approximated using the truncated Fourier series to obtain the unsteady harmonic balance equation. The implicit upper and lower symmetric Gauss-Seidel method is used to iteratively solve the equation. The inner and outer flow field variables, aerodynamic performance parameters, and aeroelastic performance parameters at different operating points are obtained, and the aerodynamic-aeroelastic multidisciplinary coupling objective function is calculated. Step 4: Based on the adjoint principle and the unsteady harmonic balance equation obtained in step 3, establish the unsteady adjoint equation and iteratively solve it to obtain the unsteady adjoint variable; Step 5: Based on the optimized design parameters obtained in step 2, the inner and outer flow field variables and the aerodynamic-aeroelastic multidisciplinary coupling objective function obtained in step 3, and the unsteady adjoint variables obtained in step 4, the partial derivative terms in the sensitivity information are calculated using the second-order central difference scheme to obtain the sensitivity information of the aerodynamic-aeroelastic multidisciplinary coupling objective function with respect to the design parameters. Step 6: Based on the sensitivity information obtained in step 5, set the perturbation step size through the line search method, and use the steepest descent method to update the optimization design parameters; Step 7: Superimpose the optimized design parameters obtained in step 6 on the blade modeling parameters in step 1, update the blade modeling parameters and blade geometry, and update the computational grid using the linear elastic method; Step 8: Replace the blade shaping parameters of step 1 with the blade shaping parameters obtained in step 7, and repeat steps 2 to 7. When the two-norm of the sensitivity information meets the convergence criterion, output the optimized blade shaping parameters and blade geometry.

[0008] Furthermore, the two-dimensional leaf profile of the leaf height extracted in step 1 needs to ensure that the leaf profiles at the leaf root, leaf midpoint, and leaf tip are all considered; usually, these seven leaf heights are 10%, 20%, 35%, 50%, 65%, 80%, and 90% respectively; The blade shaping parameters of the two-dimensional blade include the maximum deflection , relative position of maximum disturbance , relative position of maximum thickness , inlet airflow angle , outlet airflow angle , Blade installation angle and Nurbs curve control point coordinates 、 and .

[0009] Furthermore, the expression of the Hicks-Henne type function in step 2 is: Where m represents the number of design parameters corresponding to each blade shaping parameter; represents the i-th design parameter corresponding to the j-th blade shaping parameter; is the perturbation of the j-th blade modeling parameter; is the dimensionless position of the blade in the axial direction, which is expressed as follows: is the radial dimensionless coordinate, defined as follows: Where r is the radial coordinate; Indicates the radial coordinate of the hub; Represents the radial coordinate of the wheel rim.

[0010] Furthermore, the expression of the unsteady flow field equation in step three is: Where Q is the unsteady flow field variable, and R is the spatial residual of the equation; The unsteady flow field variables are approximated based on the truncated Fourier series, and the specific expression is: Where t is the physical time; is the time-averaged flow field variable; and is the i-th pair of Fourier coefficients of the flow field variables; N is the number of harmonics; is the i-th angular frequency; The specific implementation steps of establishing the unsteady harmonic balance equation are: First, solve the partial derivatives of the unsteady flow field variables with respect to time, and the expression of the time spectrum source term can be obtained as: Where, is the time partial derivative operator; is the Fourier coefficient matrix of the unsteady flow field variables; is the inverse matrix of discrete Fourier transform; E is the time spectrum source term operator; Next, substituting the time spectrum source term into the unsteady flow field equation, the unsteady harmonic balance equation can be obtained: is the spatial residual term; The specific implementation steps of using the implicit upper and lower symmetric Gauss-Seidel iteration method to solve the unsteady harmonic balance equation are as follows: First, add the pseudo-time iteration term to the left side of the unsteady harmonic balance equation, and we get: Where, represents pseudo time; Then, the spatial residual term is discretized implicitly and the time spectral source term is discretized explicitly, and the time-discrete unsteady harmonic balance equation is obtained: Where, is the pseudo-time iteration step; the superscript k represents the pseudo-time iteration step; is the flow field variable at step k; is the flow field variable at step k+1; is the time spectrum source term of the kth step; is the spatial residual term of the k+1th step; Then, the Taylor series is used to expand the spatial residual term and its linear part is taken to obtain: Where, is the increment of the unsteady flow field variable; A is the Jacobian matrix; Secondly, substitute the linearized spatial residual term into the time-discrete unsteady harmonic balance equation, move the terms and merge the similar terms, and we can get: I is the identity matrix; The coefficient matrix of the above equation is split into a diagonal matrix B, an upper triangular matrix U and a lower triangular matrix L, factorized and high-order terms are neglected, and the following equation is obtained: Next, the symmetric Gauss-Seidel iteration method is used to solve the unsteady harmonic balance equation to obtain the increments of the flow field variables.

[0011] The unsteady flow field equation solving process needs two steps of iteration, which are: First step: forward scanning; In the formula, the superscript 1 / 2 represents the forward scanning; Second step: backward scanning; Finally, the unsteady flow field variables are updated based on the increments of the unsteady flow field variables, and the following equation is obtained: The internal and external bypass aerodynamic performance parameters of different operating points include mass flow , pressure ratio , efficiency and aerodynamic performance parameter accumulated work ; subscript i represents different bypasses: b represents external bypass, and c represents internal bypass; subscript j represents different operating points: st represents near stall point, pe represents highest efficiency point, and ch represents blockage point; the calculation formulas of the above aerodynamic and aerodynamic performance parameters are as follows: In the formula, is the density; is the normal velocity perpendicular to the control surface; ds is the control unit area; is the inlet total pressure; is the outlet total pressure; is the temperature ratio; is the unsteady aerodynamic load on the blade surface; is the grid motion speed; is the unit external normal direction on the blade surface; The calculation of the aerodynamic-aerodynamic multidisciplinary coupling objective function needs two steps: First, the aerodynamic performance of the internal and external bypasses is used as the constraint, and the penalty function method is used to add the constraint, and the specific expression is as follows: Wherein, subscript 0 represents initial design; is a mass flow rate penalty function coefficient; is a total pressure ratio penalty function coefficient; is an isentropic efficiency penalty function coefficient; is an objective function under a certain working condition; is a mass flow rate of an original blade profile; is a pressure ratio of an original blade profile; is an isentropic efficiency of an original blade profile; is an accumulated work of an original blade profile; i represents different ducts: b represents outer duct, and c represents inner duct; Then, the objective functions of different working points are weighted to obtain an aerodynamic-aeroelastic multidisciplinary coupling optimization objective function in a weighted manner, and a specific expression is as follows: Wherein, is an objective function under a certain working condition; I is an objective function; The aerodynamic-aeroelastic multidisciplinary coupling objective function is expressed in a symbolic form as follows: Wherein, is a design parameter.

[0012] Further, the specific implementation steps of the step of establishing the unsteady adjoint equation are as follows: First, the unsteady harmonic balance equation of step three is linearized, and the following equation can be obtained: Secondly, the linearized objective function can be obtained as follows: Then, the linearized unsteady harmonic balance equation is substituted into the linearized objective function to obtain a sensitivity calculation formula: Finally, the is defined as an unsteady adjoint variable , and the expression of the unsteady adjoint equation is as follows: The unsteady adjoint equation is iteratively solved to obtain an unsteady adjoint variable .

[0013] Further, the sensitivity calculation formula of the aerodynamic-aeroelastic multidisciplinary coupling objective function with respect to the design parameter in the step five is as follows: The partial derivative term in the sensitivity information is calculated by using a second-order central difference format, and a specific expression is as follows: In the formula, is the disturbance quantity of the design parameter.

[0014] Further, the specific process of step six for setting the disturbance step size by the line search method is as follows: on the basis of the initial design and the objective function , an arbitrary disturbance step size is given, then the sensitivity information obtained in step five is used to update the design parameter . The update of the design parameter adopts the steepest descent method, and the specific expression is as follows: The non-steady harmonic balance equation is solved again by using the updated design parameter, the objective function is updated, and the updated objective function is compared with the initial objective function ; if the objective function value decreases, that is, , the optimization process of the next step is continued, otherwise, the disturbance step size is reduced by half, the design parameter is updated again, and the process is repeated until the objective function value decreases.

[0015] Further, the expression for updating the calculation grid by the linear elasticity method in step seven is as follows: In the formula, the subscript ss represents the suction surface; ps represents the pressure surface; dX represents the coordinate disturbance quantity of the grid point; and r represents the dimensionless distance from the suction surface to the intermediate position. Further, the sensitivity two-norm convergence criterion in step eight is as follows: .

[0016] The beneficial effects of the present application are as follows: The present application parameterizes the two-dimensional blade profile by using the maximum deflection and the relative position thereof, the maximum thickness relative position, the inlet and outlet airflow angles, the blade installation angle and the Nurbs curve control point coordinates, and fits the radial disturbance quantity of the blade modeling parameter by using the Hicks-Henne type function. Compared with the traditional blade parameterization method, the blade parameterization model proposed in the present application has more optimization degrees of freedom, a larger optimization design space and stronger practicability of the optimized blade profile. The present application calculates the sensitivity information of the objective function about the design parameter by using the non-steady companion model. The calculation model can obtain the sensitivity information of the objective function about all design parameters by solving the flow field and the companion field once respectively, and is independent of the number of design parameters. Compared with the conventional finite difference method, the present application can greatly reduce the calculation resource consumption of the sensitivity analysis and improve the optimization efficiency. The present invention adopts a penalty function method to separately constrain the inner and outer bypass aerodynamic performance parameters of the large bypass ratio fan blades, and adopts weighted objective functions of different operating points to obtain the overall objective function. Compared with the existing multidisciplinary coupling optimization technology of single-duct fan blades, it can improve the aerodynamic performance of the large bypass ratio fan in the full operating range while ensuring that the aerodynamic performance of the inner and outer bypass at different operating points is not reduced, providing important support for the subsequent development of multidisciplinary coupling optimization technology for large bypass ratio fan blades. BRIEF DESCRIPTION OF THE DRAWINGS

[0017] Figure 1 This is a schematic diagram of the multidisciplinary coupling optimization process for high bypass ratio fan blades of the present invention.

[0018] Figure 2 Schematic diagram of the meridian plane of the computational mesh for a high bypass ratio fan.

[0019] Figure 3 Schematic diagram of the distribution of leaf shaping parameters at different leaf heights.

[0020] Figure 4 Schematic diagram of three-dimensional Hicks-Henne function distribution.

[0021] Figure 5 Schematic diagram of the convergence history for the optimization of aerodynamic and aeroelastic performance parameters.

[0022] Figure 6 Schematic diagram for comparison of aeroelastic performance parameters before and after optimization.

[0023] Figure 7 Schematic diagram comparing the aerodynamic performance parameters of the inner and outer ducts before and after optimization.

[0024] Figure 8 Schematic diagram for comparing leaf profiles with different leaf heights before and after optimization. DETAILED DESCRIPTION

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

[0026] In order to achieve the above object, the present invention provides the following specific implementation methods: Example 1: Figure 1 As shown, a multidisciplinary coupled adjoint optimization method for a high bypass ratio fan blade includes the following steps: Step 1: Extract the two-dimensional blade profiles of seven different blade heights, 10%, 20%, 35%, 50%, 65%, 80%, and 90%, along the span direction of the high bypass ratio fan blade; for each blade height, extract the blade modeling parameters, including the maximum deflection t, the relative position of the maximum deflection , relative position of maximum thickness , inlet airflow angle outlet flow angle blade mounting angle and Nurbs curve control point coordinates , , and ; through this step, two-dimensional blade modeling parameters of different blade heights can be extracted and applied to the parameterization of the radial disturbance amount of the modeling parameter in the next step; Step two: the radial disturbance amount of the blade modeling parameter obtained in step one is parameterized by using the Hicks-Henne type function method to obtain the optimization design parameter, and the specific expression of the Hicks-Henne type function is as follows: In the formula, m represents the number of design parameters corresponding to each blade modeling parameter; subscript j represents different blade modeling parameters; represents the i-th design parameter corresponding to the j-th blade modeling parameter; is the disturbance amount of the j-th blade modeling parameter; is defined as follows: is a radial dimensionless coordinate, and is defined as follows: Wherein, r is a radial coordinate; subscript hub represents a hub; and subscript shroud represents a shroud; the parameterization of the radial disturbance amount of the blade modeling parameter by the Hicks-Henne type function can ensure that the optimized blade profile is smooth along the spanwise direction; Step three: based on the unsteady flow field equation, the unsteady harmonic balance equation is obtained by approximating the unsteady flow field variable based on the truncated Fourier series, and the implicit upper and lower symmetric Gauss-Seidel method is used for iteration to obtain the internal and external flow field variables, aerodynamic performance parameters and aeroelastic performance parameters of different working points, and to calculate the aerodynamic-aeroelastic multidisciplinary coupling objective function. The specific expression of the unsteady flow field equation is as follows: In the formula, Q is an unsteady flow field variable; and R is an equation space residual; The specific implementation steps of establishing the unsteady harmonic balance equation are as follows: First, the unsteady flow field variable is approximated based on the truncated Fourier series, and the specific expression is as follows: Secondly, the partial derivative of the unsteady flow field variable with respect to time is solved, and the time spectrum source term expression is as follows: In the formula, is the time derivative operator; is the Fourier coefficient matrix of the unsteady flow field variables; is the inverse discrete Fourier transform matrix; E is the time spectrum source term operator; Finally, substituting the time spectrum source term into the unsteady flow field equation, the unsteady harmonic balance equation is obtained: The specific implementation steps of solving the unsteady harmonic balance equation by using the implicit upper and lower symmetric Gauss-Seidel iteration method are as follows: First, a pseudo-time iteration term is added to the left side of the unsteady harmonic balance equation, and the following equation is obtained: In the formula, represents the pseudo-time; Next, the spatial residual term is discretized by using the implicit method and the time spectrum source term is discretized by using the explicit method, and the time-discrete unsteady harmonic balance equation is obtained: In the formula, is the pseudo-time iteration step; the superscript k represents the pseudo-time iteration step; Then, the spatial residual term is expanded by using the Taylor series and the linear part is taken, and the following equation is obtained: In the formula, is the increment of the unsteady flow field variables; A is the Jacobian matrix; Secondly, the linearized spatial residual term is substituted into the time-discrete unsteady harmonic balance equation, and the same terms are removed and combined, and the following equation is obtained: The coefficient matrix of the above equation is split into a diagonal matrix B, an upper triangular matrix U and a lower triangular matrix L, factorized and ignored high-order terms, and the following equation is obtained: Next, the upper and lower symmetric Gauss-Seidel iteration method is used to solve the unsteady harmonic balance equation to obtain the increment of the flow field variables. The unsteady flow field equation solving process needs two steps of iteration, which are: First step: forward scanning; In the formula, the superscript 1 / 2 represents the forward scanning; Second step: backward scanning; Finally, based on the increment of the unsteady flow field variables, the unsteady flow field variables are updated, and the following equation is obtained: Based on the above unsteady flow field variables, the internal and external duct aerodynamic performance parameters including mass flow , pressure ratio , efficiency and the aerodynamic and aeroelastic performance parameters accumulated work of different operating points are calculated. Subscript i represents different ducts: b represents the outer duct, c represents the inner duct; subscript j represents different operating points: st represents the near stall point, pe represents the highest efficiency point, ch represents the blockage point. The calculation formulas of the above aerodynamic and aeroelastic performance parameters are as follows: In the formula, is the density; is the normal velocity perpendicular to the control surface; ds is the control unit area; is the inlet total pressure; is the outlet total pressure; is the temperature ratio; is the unsteady aerodynamic load on the blade surface; is the grid motion velocity; is the unit outer normal direction on the blade surface; According to the internal and external duct aerodynamic performance parameters and the aeroelastic performance parameters of different operating points, the aerodynamic-aeroelastic multidisciplinary coupling objective function is calculated, and the specific implementation process is as follows: First, taking the aeroelastic performance parameters as the target and the internal and external duct aerodynamic performance as the constraint, the penalty function method is used to add the constraint, and the specific expression is as follows: In the formula, subscript 0 represents the initial design; is the mass flow penalty function coefficient; is the total pressure ratio penalty function coefficient; is the isentropic efficiency penalty function coefficient; is the objective function under a certain operating point; Then, the objective functions of different operating points are weighted and processed in a weighted manner to obtain the aerodynamic-aeroelastic multidisciplinary coupling optimization objective function, and the specific expression is as follows: In the formula, is the weighting coefficient of the objective function under a certain operating point; I is the objective function; The above aerodynamic-aeroelastic multidisciplinary coupling objective function can be expressed in the following symbolic form: In the formula, is the design parameter; through this step, the unsteady flow field variables and the objective function can be calculated, which are used to solve the adjoint variables in the next step; Step four: according to the unsteady harmonic balance equation obtained in step three and the adjoint principle, the unsteady adjoint equation is established, and the unsteady adjoint variables are obtained by iteration. The specific implementation steps of establishing the unsteady adjoint equation are as follows: Firstly, the unsteady harmonic balance equation in step three is linearized, and the following equation can be obtained: Secondly, the linearized objective function can be obtained as follows: Then, the linearized unsteady harmonic balance equation is substituted into the linearized objective function to obtain the sensitivity calculation formula: Finally, let be the unsteady adjoint variable , then the expression of the unsteady adjoint equation is as follows: The unsteady adjoint equation is solved by iteration, and the unsteady adjoint variable is obtained, which is used to solve the sensitivity information of the objective function with respect to the design parameter in the next step; Step five: based on the optimization design parameter obtained in step two, the internal and external flow field variables obtained in step three, the aerodynamic-aeroelastic multi-disciplinary coupling objective function and the unsteady adjoint variable obtained in step four, the partial derivative term in the sensitivity information is calculated by using the second-order central difference format, and the sensitivity of the aerodynamic-aeroelastic multi-disciplinary coupling objective function with respect to the design parameter is obtained. The sensitivity calculation formula of the aerodynamic-aeroelastic multi-disciplinary coupling objective function with respect to the design parameter is as follows: The calculation of and in the above aerodynamic-aeroelastic multi-disciplinary coupling objective function adopts the second-order central difference format, and the specific calculation formula is as follows: In the formula, is the perturbation of the design parameter; through this step, the sensitivity information of the objective function with respect to the design parameter is obtained, which is used to update the design parameter and the blade profile in the next step; Step six: based on the sensitivity information obtained in step five, the perturbation step is set by the line search method, and the optimization design parameter is updated by using the steepest descent method. The specific implementation process of the line search method is as follows: in the initial design and the objective function on the basis of the given perturbation step and then the sensitivity information obtained in step five , the design parameters are updated The update of the design parameters uses the steepest descent method, and the specific expression is as follows: Based on the updated design parameters, the unsteady harmonic balance equation is solved again, and the objective function is updated, which is compared with the initial objective function If the objective function value decreases, that is , the next optimization process is continued, otherwise, the perturbation step is reduced by half, the design parameters are updated again, and the optimization task of updating the design parameters and realizing the decrease of the objective function value is completed through this step; Step seven: The optimized design parameters obtained in step six are superimposed on the blade modeling parameters in step one, the blade modeling parameters and the blade geometry are updated, and the linear elastic method is used to update the calculation grid. The specific expression of the linear elastic method is: In the formula, the subscript ss represents the suction surface; ps represents the pressure surface; dX represents the coordinate perturbation of the grid point; r represents the dimensionless distance from the suction surface to the middle position; for small deformation optimization problems, the grid perturbation method based on linear elasticity not only ensures the quality of the optimized grid, but also is more efficient; Step eight: The blade modeling parameters obtained in step seven are replaced with the blade modeling parameters in step one, and steps two to seven are repeated, and when the two-norm of the sensitivity meets the convergence criterion, the optimized blade modeling parameters and blade geometry are output, and the specific expression of the two-norm convergence criterion of the sensitivity is: .

[0027] Through this step, the optimization termination condition can be determined without human intervention, and the automation of the entire optimization process is realized.

[0028] In order to more specifically describe the purpose, technical scheme and advantages of the present application, a large-bypass-ratio fan is taken as an example for further detailed description, Figure 2 The meridian plane schematic diagram of the large-bypass-ratio fan is shown, including a test inlet, a test inner bypass outlet and a test outer bypass outlet, and the mass flow ratio of the inner and outer bypasses is 1:9: Example: As Figures 2 to 8As shown, the flow field solver and the accompanying solver of the present example are based on self-programmed Fortran programs, the optimization platform is based on self-programmed python script programs, the blade modeling adopts self-programmed T-blade programs, the radial disturbance parameterization of the two-dimensional blade profile parameters adopts self-programmed Fortran programs, and the grid disturbance based on linear elasticity is based on self-programmed Fortran programs. The specific implementation steps of the present example are as follows: Step (1): Based on the self-programmed blade modeling program T-blade, first extract the two-dimensional blade profile at seven different blade heights of 10%, 20%, 35%, 50%, 65%, 80% and 90%. For the two-dimensional blade profile at each blade height, extract its ten blade modeling parameters, including the maximum deflection , the maximum disturbance relative position , the maximum thickness relative position , the inlet flow angle , the outlet flow angle , the blade installation angle and the Nurbs curve control point coordinates , and , and write the blade modeling parameters into the row.prof file, Figure 3 The distribution diagram of the two-dimensional blade modeling parameters at the blade height of 50% is shown in FIG. 1; Step (2): Based on the self-programmed blade parameterization program, read the row.prof file in step (1), for each two-dimensional blade modeling parameter, parameterize the radial disturbance amount by using ten Hicks-Henne type functions, write the design parameter information of the Hicks-Henne type function, including the number of modeling parameters and the number of design parameters 100, into the dgn_var_inf.dat file, write the initial value information of each design parameter into the dgn_var.dat file, and the initial value of the design parameter is usually set to 0; as shown in FIG. 2; Figure 4 ​Step (3): Fluid mesh is divided based on Autogrid5 module in NUMECA software, and.grd mesh file is saved, and blade vibration mode and vibration frequency information is calculated based on Workbench module in ANSYS software, and.in mode file is saved. Based on self-compiled flow field solver,.grd mesh file and.in mode file are read respectively, fluid-structure coupling analysis is carried out, unsteady flow field variables are saved in solution.save file, and aerodynamic performance parameters including mass flow, pressure ratio and efficiency of inner and outer ducts and aeroelastic performance parameters including accumulated work of three different working conditions including near stall, blockage and highest efficiency are calculated. The initial aerodynamic performance parameters and aeroelastic performance parameters are written into objfun0.dat file. With accumulated work as objective function, mass flow, pressure ratio and efficiency of inner and outer ducts as constraints, aerodynamic-aeroelastic multi-disciplinary coupling objective function is calculated. Among them, the objective function weighting coefficients of near stall condition, blockage condition and highest efficiency point are set to 0.5, 0.25 and 0.25 respectively, and the penalty function coefficients of mass flow, pressure ratio and efficiency of inner and outer ducts are all set to 200; Step (4): Based on the self-compiled adjoint solver program, solution.save file in step (3) is read for initializing unsteady flow field variables, solving unsteady adjoint equations, obtaining unsteady adjoint variables, and saving the results in solution_adj.save; Step (5): Based on the optimization design parameters in step (2), the inner and outer duct flow field variables in step (3), the aerodynamic-aeroelastic multi-objective coupling objective function and the unsteady adjoint variables in step (4), the partial derivative terms in the sensitivity information are calculated by using the second-order central difference format, the sensitivity of the aerodynamic-aeroelastic multi-disciplinary coupling objective function with respect to the design parameters is obtained, and the results are saved in sensitivity.dat file; Step (6): By using the self-compiled python optimization script program, the sensitivity information obtained in step (5) is read, the perturbation step is set by using line search method, and the optimization design parameters are updated by using steepest descent method, and the design parameter information in dgn_var.dat is updated; Step (7): The optimization design parameters obtained in step six are superimposed on the blade modeling parameters in step one, the blade modeling parameters and blade geometry are updated, the calculation mesh is updated by using linear elastic method, and the mesh information is written into.grd file; Step (8): The blade modeling parameters in step (7) are replaced based on the blade modeling parameters in step (1), and steps (2) to (7) are repeated, when the two norm of sensitivity meets the convergence criterion, that is , the optimized blade modeling parameters and blade geometry are output.

[0029] Figure 5The evolution of the aeroelastic performance objective function and the aerodynamic performance constraints with the iteration step number in the optimization process is shown. It can be seen that after 10 optimization steps, the aeroelastic performance objective function decreases by about 8%, while the aerodynamic performance parameter constraints do not change substantially. Figure 6 The aeroelastic performance objective functions at different pitch diameters before and after optimization are compared. After optimization, the minimum logarithmic decay rate increases, and the aeroelastic stability is improved. Figure 7 The internal and external ratio of pressure and efficiency before and after optimization are compared. It can be seen that after optimization, the aerodynamic performance does not decrease. Figure 8 The blade profiles at different blade heights before and after optimization are compared. It can be seen that the reduction of blade curvature is beneficial to improve the aeroelastic stability.

[0030] The above only describes the preferred embodiments of the present application, and is not intended to limit the present application. Any modification, equivalent replacement, improvement, etc. made within the spirit and principles of the present application shall be included in the protection scope of the present application.

Claims

1. A multidisciplinary coupled adjoint optimization method for high bypass ratio fan blades, characterized in that: The following steps are included: Step 1: Extract the two-dimensional blade profiles of at least seven different blade heights of the high bypass ratio fan blade, and extract the blade shaping parameters of each blade height; Step 2: Use the Hicks-Henne function method to parameterize the radial perturbation of the blade modeling parameters obtained in step 1 to obtain the optimized design parameters; Step 3: Based on the unsteady flow field equation, the unsteady flow field variables are approximated using the truncated Fourier series to obtain the unsteady harmonic balance equation. The implicit upper and lower symmetric Gauss-Seidel method is used to iteratively solve the equation. The inner and outer flow field variables, aerodynamic performance parameters, and aeroelastic performance parameters at different operating points are obtained, and the aerodynamic-aeroelastic multidisciplinary coupling objective function is calculated. Step 4: Based on the adjoint principle and the unsteady harmonic balance equation obtained in step 3, establish the unsteady adjoint equation and iteratively solve it to obtain the unsteady adjoint variable; Step 5: Based on the optimized design parameters obtained in step 2, the inner and outer flow field variables and the aerodynamic-aeroelastic multidisciplinary coupling objective function obtained in step 3, and the unsteady adjoint variables obtained in step 4, the partial derivative terms in the sensitivity information are calculated using the second-order central difference scheme to obtain the sensitivity information of the aerodynamic-aeroelastic multidisciplinary coupling objective function with respect to the design parameters. Step 6: Based on the sensitivity information obtained in step 5, set the perturbation step size through the line search method, and use the steepest descent method to update the optimization design parameters; Step 7: Superimpose the optimized design parameters obtained in step 6 on the blade modeling parameters in step 1, update the blade modeling parameters and blade geometry, and update the computational grid using the linear elastic method; Step 8: Replace the blade shaping parameters of step 1 with the blade shaping parameters obtained in step 7, and repeat steps 2 to 7. When the two-norm of the sensitivity information meets the convergence criterion, output the optimized blade shaping parameters and blade geometry.

2. The multidisciplinary coupled adjoint optimization method for high bypass ratio fan blades according to claim 1, characterized in that: The two-dimensional leaf profiles of different leaf heights extracted in step 1 need to ensure that the leaf profiles at the root, middle, and tip are all considered; The blade shaping parameters of the two-dimensional blade include the maximum deflection , relative position of maximum disturbance , relative position of maximum thickness , inlet airflow angle , outlet airflow angle , Blade installation angle and Nurbs curve control point coordinates 、 and .

3. The multidisciplinary coupled adjoint optimization method for high bypass ratio fan blades according to claim 2, characterized in that: The expression of the Hicks-Henne type function in step 2 is: Where m represents the number of design parameters corresponding to each blade shaping parameter; represents the i-th design parameter corresponding to the j-th blade shaping parameter; is the disturbance of the j-th blade modeling parameter; is the dimensionless position of the blade in the axial direction, which is expressed as follows: is the radial dimensionless coordinate, defined as follows: Where r is the radial coordinate; Indicates the radial coordinate of the hub; Represents the radial coordinate of the wheel rim.

4. The multidisciplinary coupled adjoint optimization method for high bypass ratio fan blades according to claim 3 is characterized in that: The expression of the unsteady flow field equation in step 3 is: Where Q is the unsteady flow field variable; R is the spatial residual of the equation; The unsteady flow field variables are approximated based on the truncated Fourier series. The specific expression is: Where t is the physical time; is the time-averaged flow field variable; and is the i-th pair of Fourier coefficients of the flow field variables; N is the number of harmonics; is the i-th angular frequency; The specific implementation steps for establishing the unsteady harmonic balance equation are: First, solve the partial derivatives of the unsteady flow field variables with respect to time, and the expression of the time spectrum source term is obtained as: Where, is the time partial derivative operator; is the Fourier coefficient matrix of the unsteady flow field variables; is the inverse matrix of discrete Fourier transform; E is the time spectrum source term operator; Next, the time spectrum source term is substituted into the unsteady flow field equation to obtain the unsteady harmonic balance equation: is the spatial residual term; The specific implementation steps of the implicit upper and lower symmetric Gauss-Seidel iteration method to solve the unsteady harmonic balance equation are as follows: First, add the pseudo-time iteration term to the left side of the unsteady harmonic balance equation, and we get: Where, represents pseudo time; Then, the implicit method is used to discretize the spatial residual term and the explicit method is used to discretize the time spectral source term, and the time-discrete unsteady harmonic balance equation is obtained: Where, is the pseudo-time iteration step; the superscript k represents the pseudo-time iteration step; is the flow field variable at step k; is the flow field variable at step k+1; is the time spectrum source term of the kth step; is the spatial residual term of the k+1th step; Then, the Taylor series is used to expand the spatial residual term and its linear part is taken to obtain: Where, is the increment of the unsteady flow field variable; A is the Jacobian matrix; Secondly, substitute the linearized spatial residual term into the time-discrete unsteady harmonic balance equation, move the terms and merge the similar terms, and we get: I is the identity matrix; Split the coefficient matrix of the above equation into a diagonal matrix B, an upper triangular matrix U, and a lower triangular matrix L, factorize them and ignore the high-order terms, and we get: Next, the upper and lower symmetric Gauss-Seidel iteration method is used to solve the unsteady harmonic balance equation to obtain the increment of the flow field variables.

5. The multidisciplinary coupled adjoint optimization method for high bypass ratio fan blades according to claim 4, characterized in that: The solution process of the unsteady flow field equation requires two steps of iteration, namely: Step 1: Scan forward; Where, the superscript 1 / 2 indicates forward scanning; Step 2: Scan backward; Finally, based on the increment of the above unsteady flow field variables, the unsteady flow field variables are updated, and we get: The aerodynamic performance parameters of the inner and outer ducts at different working points include mass flow rate ; Pressure ratio ;efficiency and aeroelastic performance parameter accumulation work The subscript i represents different ducts: b represents external ducts, and c represents internal ducts. The subscript j represents different operating points: st represents the near-stall point, pe represents the highest efficiency point, and ch represents the clogging point. The calculation formulas for the above aerodynamic and aeroelastic performance parameters are: Where, is the density; is the normal velocity perpendicular to the control surface; ds is the control unit area; is the total inlet pressure; is the total outlet pressure; For temperature ratio; is the unsteady aerodynamic load on the blade surface; is the grid movement speed; is the unit normal direction of the blade surface.

6. The multidisciplinary coupled adjoint optimization method for high bypass ratio fan blades according to claim 5, characterized in that: The calculation of the aerodynamic-aeroelastic multidisciplinary coupling objective function requires two steps: First, the aeroelastic performance parameters are used as the target, and the internal and external aerodynamic performance are used as the constraints. The constraint is added using the penalty function method. The specific expression is as follows; Where, subscript 0 indicates the initial design; is the mass flow penalty function coefficient; is the total pressure ratio penalty function coefficient; is the coefficient of the isentropic efficiency penalty function; is the objective function under a certain working condition; is the mass flow rate of the original blade; is the pressure ratio of the original blade profile; is the isentropic efficiency of the original blade type; is the accumulated work of the original blade type; i represents different ducts: b represents external duct, c represents internal duct; Then, the objective functions of different working points are weighted in a weighted manner to obtain the aerodynamic-aeroelastic multidisciplinary coupling optimization objective function. The specific expression is as follows: Where, is the weighted coefficient of the objective function under a certain working condition; I is the objective function; The above-mentioned aerodynamic-aeroelastic multidisciplinary coupling objective function is expressed in the following symbolic form: Where, is the design parameter.

7. The multidisciplinary coupled adjoint optimization method for high bypass ratio fan blades according to claim 6, characterized in that: The specific implementation steps of step 4 to establish the unsteady adjoint equation are: First, linearize the unsteady harmonic balance equation in step 3 and obtain: Secondly, the linearized objective function can be obtained: Then, the linearized unsteady harmonic balance equation is substituted into the linearized objective function to obtain the sensitivity calculation formula: Finally, Defined as an unsteady adjoint variable , the expression of the unsteady adjoint equation is as follows: Iteratively solve the above unsteady adjoint equation to obtain the unsteady adjoint variable .

8. The multidisciplinary coupled adjoint optimization method for high bypass ratio fan blades according to claim 7, characterized in that: The calculation formula for the sensitivity of the aerodynamic-aeroelastic multidisciplinary coupling objective function to the design parameters in step 5 is: The second-order central difference format is used to calculate the partial derivative terms in the sensitivity information. The specific expression is: Where, is the disturbance of the design parameters.

9. The multidisciplinary coupled adjoint optimization method for high bypass ratio fan blades according to claim 8, characterized in that: Step 6 The specific process of setting the perturbation step size by line search method is as follows: and the objective function Based on the perturbation step size, , and then based on the sensitivity information obtained in step five , update the design parameters , the design parameters are updated using the steepest descent method, and the specific expression is as follows: Resolve the unsteady harmonic balance equation based on the updated design parameters and update the objective function , and compare it with the initial objective function For comparison; if the objective function value decreases, that is , then continue the optimization process to the next step, otherwise, reduce the perturbation step size by half and re-update the design parameters until the objective function value decreases.

10. The multidisciplinary coupled adjoint optimization method for high bypass ratio fan blades according to claim 9, characterized in that: Step 7: Update the computational grid using the linear elastic method. The expression is: Where, the subscript ss represents the suction surface; ps represents the pressure surface; dX represents the coordinate perturbation of the grid point; r represents the dimensionless distance from the suction surface to the middle position; The sensitivity two-norm convergence criterion of step eight is specifically expressed as follows: 。

Images