Unsteady aerodynamic adjoint optimization method for aeroaxial turbomachinery considering multi-source disturbances

Abstract

The application discloses a method for considering multi-source disturbance of an aviation axial flow impeller machine unsteady aerodynamic adjoint optimization, which comprises the following steps: step one, grid division of aviation axial flow impeller blades; step two, extraction of flow field information; step three, determination of the number of harmonics required for analyzing different disturbances; step four, calculation of unsteady flow frequency and inter-blade phase angle of different disturbances; step five, non-uniform time point sampling by using approximate periodic Fourier transformation; step six, acquisition of unsteady flow field information; step seven, solution of an adjoint equation to obtain sensitivity information of a target function about design parameters; step eight, update of the aviation axial flow impeller blade shape by using the steepest descent method; step nine, repetition of step one, step six, step seven and step eight until the optimization result meets the convergence standard, and output of the optimization result. The application has the outstanding characteristics of low calculation cost and optimization result closer to the engineering practical application requirements.

Description

Technical Field

[0001] This invention relates to the field of unsteady aerodynamic flow field analysis and blade shape optimization for aerospace axial flow turbomachinery, specifically to an unsteady aerodynamic adjoint optimization method for aerospace axial flow turbomachinery considering multi-source disturbances. Background Technology

[0002] Axial-flow turbomachinery (including fans, compressors, and turbines) is a core component of aero-engines, and its aerodynamic performance directly affects key engine parameters such as thrust-to-weight ratio and fuel consumption. With the continuous pursuit of higher aerodynamic loads and superior performance in advanced aero-engines, the internal flow of aero-axial-flow turbomachinery exhibits increasingly complex characteristics, such as shock waves, shock wave-boundary layer interference, and inter-row interference. These typical flow patterns have a significant impact on the aerodynamic performance of aero-axial-flow turbomachinery and are extremely sensitive to subtle changes in blade geometry. Against this backdrop, employing refined optimization design methods to optimize the three-dimensional blade shape of aero-axial-flow turbomachinery has become an important approach to reduce flow losses and further unlock the aerodynamic performance potential of aero-axial-flow turbomachinery.

[0003] In aerodynamic optimization design methods, the adjoint method based on control theory has attracted widespread attention due to its unique advantage that the computational cost of solving the gradient of the objective function with respect to design variables is essentially independent of the number of design variables. Compared with traditional gradient solving methods, the adjoint method can implement refined three-dimensional blade design with lower computational cost even with numerous design parameters, making it an important development direction for aerodynamic design technology of aerospace axial turbines. Currently, scholars both domestically and internationally have conducted extensive research on the adjoint method and achieved significant results in the aerodynamic optimization of typical turbine components such as axial turbine blades and fan / compressor blades.

[0004] However, existing adjoint optimization methods suffer from two main shortcomings: one type is based on the assumption of a steady flow field, failing to reflect the unsteady nature of actual flow, as exemplified by the paper titled "Multi-Row Turbomachinery Aerodynamic Design Optimization by an Efficient and Accurate Discrete Adjoint Method"; the other type, while introducing unsteady effects, only analyzes and optimizes a single disturbance source, as exemplified by the paper titled "Three-Dimensional Multi-Row Unsteady Aerodynamic Optimization of a Compressor Stage Using an Implicit Discrete Adjoint Harmonic Balance Method". In reality, aero-engines inevitably experience multiple disturbance sources simultaneously during actual operation. These disturbances include, but are not limited to, inlet distortion, interference between dynamic and static exhaust stages, manufacturing tolerances, and changes in the operating environment. Multi-source disturbances induce complex unsteady responses in the internal flow field of the turbine, and the unsteady effects caused by different disturbance sources couple and superimpose, resulting in the actual aerodynamic performance of the blades deviating significantly from the design expectations. Currently, there are no publicly available technical disclosures regarding turbomachinery adjoint optimization methods that simultaneously consider the aforementioned multi-source disturbances.

[0005] In summary, existing aerodynamic optimization techniques for aero-axial flow turbofans fail to effectively consider the multi-source disturbances present in actual operation, making it difficult to meet the high reliability and robustness requirements of advanced aero-engine aerodynamic optimization design. Therefore, there is an urgent need to develop an aerodynamic optimization method for aero-axial flow turbofans that can simultaneously consider multi-source disturbances and unsteady flow characteristics, in order to overcome the shortcomings of existing technologies and provide more reliable technical support for the refined and robust aerodynamic design of aero-axial flow turbofans. Summary of the Invention

[0006] To overcome the shortcomings of existing technologies, this invention proposes an unsteady aerodynamic adjoint optimization method for aerospace axial-flow turbofans that considers multi-source disturbances. This method simultaneously considers two different types of aerodynamic disturbances: upstream wake and downstream potential disturbances, enabling highly reliable and robust optimization of turbofan aerodynamic performance within an unsteady framework. Furthermore, this invention fully utilizes the gradient-efficient solution characteristics of the adjoint method, resulting in low computational resource consumption and optimization results that more closely reflect the needs of practical engineering applications.

[0007] To achieve the above objectives, the technical solution adopted by the present invention is as follows: An unsteady aerodynamic adjoint optimization method for aerospace axial flow turbines considering multi-source disturbances includes the following steps; Step 1: Perform fluid mesh generation on the internal flow channels of the aerospace axial flow turbine, and structure the topology of the mesh generation into a single block to obtain the fluid computation mesh required for flow field analysis; Step 2: Based on the fluid computational grid obtained in Step 1, solve the steady Navier-Stokes equations and extract the turbulence variable cloud map and its circumferential distribution at 50% leaf height; Step 3: Based on the circumferential distribution of turbulent variables at 50% blade height in Step 2, perform a fast Fourier transform to determine the number of harmonics required to resolve different disturbances in the unsteady flow field analysis. Step 4: Based on the design speed, operating conditions, and number of blades in different blade rows of the aero-axial flow turbine shaft, determine the unsteady flow frequency of different disturbances and the inter-blade phase angle information corresponding to each frequency component; Step 5: Based on the number of harmonics obtained in Step 3 and the unsteady current frequency obtained in Step 4, non-uniform time point sampling is performed using the approximate periodic Fourier transform technique; Step Six: Based on the semi-discrete form of the unsteady Navier-Stokes equations, and using the number of harmonics obtained in Step Three, the unsteady flow frequency obtained in Step Four, the non-uniform time points obtained in Step Five, and the principle of harmonic balance technology, construct the unsteady harmonic balance control equations, and solve them discretely to obtain the unsteady flow field information inside the aerospace axial flow turbine. Step 7: Input the unsteady flow field information obtained in Step 6 into the adjoint optimization solver, solve the adjoint equation, obtain the adjoint sensitivity information of the objective function with respect to the design parameters, and verify the accuracy of the adjoint sensitivity calculation using the finite difference method. Step 8: Update the shape of the aerospace axial flow impeller blades based on the sensitivity information of the objective function with respect to design parameters and the steepest descent method described in Step 7; Step 9: Repeat steps 1, 6, 7 and 8 until the optimization result meets the convergence criteria.

[0008] Furthermore, in step one, the impeller blade fluid mesh generation uses the Autogrid5 module of the commercial software NUMECA, and the mesh topology is HI monolithic structured.

[0009] Furthermore, the expression for the steady Navier-Stokes equations in step two is as follows: In the formula, For flow field variables; Density; U is the velocity vector; For internal energy; For turbulence variables; F, G, and H are the convective fluxes in the axial, circumferential, and radial directions, respectively. , and , respectively, represent the viscous fluxes in the axial, circumferential, and radial directions; S is the source term; For pseudo-time; x, r and r are the axial, circumferential, and radial coordinates, respectively.

[0010] Solving the steady Navier-Stokes equations in step two yields the flow field variable Q within the entire fluid grid, and Q contains... For turbulent variables, extract the flow field variables at 50% of the blade height within the entire fluid grid to obtain the turbulent variable contour map at 50% of the blade height. Based on the turbulent variable contour map at 50% of the blade height, extract the turbulent variable distribution at the outlet to obtain the circumferential distribution of the turbulent variables at 50% of the blade height.

[0011] Furthermore, the specific implementation process of the Fast Fourier Transform described in step three is as follows: First, the circumferential distribution of the turbulent variables at 50% of the blade height obtained in step two is expressed as: In the formula, N is the number of circumferential grid points; The turbulence variable for the first grid point in the circumferential direction; For the turbulence variables at the second grid point in the circumferential direction; Let N be the turbulence variable at the Nth grid point in the circumferential direction; Secondly, regarding Performing a discrete Fourier transform yields: In the formula, The rotation factor; This represents the k-th Fourier coefficient; Let n be the turbulence variable at the nth circumferential grid point; Next, based on the periodicity of the Discrete Fourier Transform, the formula for the butterfly operation unit of the Fast Fourier Transform is obtained, as shown in the following expression: In the formula, the subscript j represents the order of the Fourier coefficients; N represents the number of sampling points; It is the rotation factor; Let be the j-th order Fourier coefficients; The discrete Fourier transform of N / 2 points with an even number of points; The discrete Fourier transform of N / 2 points with an odd number of points is expressed as follows: Finally, in order to determine the number of harmonics required to resolve different disturbances, the following dimensionless parameters are defined: Make If this holds true, the smallest k corresponding to this condition is the number of harmonics required to resolve different disturbances.

[0012] Furthermore, the expressions for the unsteady current frequencies of the different disturbances mentioned in step four are as follows: In the formula, The frequency is unsteady. is the design rotational speed of the shaft; i is the harmonic number; Operating conditions; Pi is the mathematical constant of a circle.

[0013] Furthermore, the expression for the inter-leaf phase angle corresponding to each frequency component described in step four is as follows: In the formula, and This indicates the number of blades in the first and second rows.

[0014] Furthermore, the specific implementation steps of the approximate periodic Fourier transform technique described in step five for sampling at non-uniform time points are as follows: First, based on the unsteady current frequency calculated in step four, the minimum unsteady current frequency information is obtained, as shown in the following formula: Secondly, based on the minimum unsteady current frequency information, the time information is divided into M equal parts, resulting in M ​​time-sampling points. The specific calculation formula is as follows: In the formula, Let M be the sampling point at time m. In actual calculations, M is usually set to 40000. Finally, based on the Gram-Schmidt orthogonalization method, 2k+1 non-uniform sampling points are obtained. The specific implementation process of the Gram-Schmidt orthogonalization method is as follows: 1) Based on the time point information, construct the column vector as follows: 2) Iterative solution minimizes the condition number of the discrete Fourier transform matrix. The specific formula is as follows: In the formula, the superscript T denotes the transpose of the vector; For time points The corresponding column vector; for The column vector corresponding to the given time; 3) Obtain 2k+1 non-uniform sampling time points, as follows: In the formula, The non-uniform sampling time points are the solution to be required.

[0015] Furthermore, the specific expression for the semi-discrete form of the unsteady Navier-Stokes equations described in step six is ​​as follows: In the formula, Q represents the flow field variable, and R represents the equation space residual; Next, the flow field variables are approximated using a truncated Fourier series, expressed as: In the formula, These are time-averaged flow field variables; and These are the Fourier coefficients of the flow field variables; The interleaf phase angle; The frequency is unsteady. For time; Pi; Then, by solving the time partial derivative of the above Fourier expression, we obtain the expression for the time spectrum source term, as follows: In the formula, For time partial derivative operators; is the Fourier coefficient vector of the flow field variables; The matrix is ​​the inverse of the discrete Fourier transform; E is the time spectrum source term operator. Finally, substituting the time spectrum source term into the unsteady flow field control equation, we obtain the aforementioned unsteady harmonic balance control equation: Solving the above unsteady harmonic equilibrium control equations yields the flow field variable Q, thus providing the unsteady flow field information.

[0016] Furthermore, the expression for the adjoint equation described in step seven is as follows: In the formula, As a companion variable; I is the transpose of the time spectrum source term operator; I is the objective function, which is the mass flow rate, total pressure ratio, isentropic efficiency, or a combination thereof. It is pseudo-time; Solve the above adjoint equation to obtain the adjoint variables. Substitute the adjoint variables into the sensitivity calculation formula to obtain the adjoint sensitivity described in step seven. The specific calculation formula is as follows: In the formula, To provide design parameters, the blade profile is parameterized using Hicks-Henne type functions. The expression for the Hicks-Henne type function is as follows: In the formula, This represents the perturbation amount of the circumferential grid coordinates; , , and The definition is as follows: In the formula, the subscripts le and te represent the leading edge and trailing edge, respectively; the subscripts hub and shroud represent the rim and hub, respectively; and M represents the radially extracted blade height number.

[0017] Furthermore, the expression for the objective function described in step seven is: In the formula, the subscript 0 represents the initial value; This is due to entropy increase; For mass flow rate; This is the mass flow penalty function coefficient, with a value of 200; Total pressure ratio; This is the total pressure ratio penalty function coefficient, with a value of 200.

[0018] Furthermore, the sensitivity calculation formula for the finite difference method described in step seven is as follows: In the formula, This represents the disturbance amount of the design parameters.

[0019] Furthermore, the specific expression for updating the design parameters using the steepest descent method described in step eight is as follows: In the formula, n represents the nth optimization step; This represents the disturbance step size, which is usually set to 1% of the blade chord length.

[0020] Furthermore, step nine optimizes the convergence criterion by setting the sensitivity's L2 norm to be less than 1. The specific expression is as follows: .

[0021] The beneficial effects of this invention are: This invention employs the adjoint method to calculate the sensitivity of the objective function of an aerospace axial-flow turbomachinery to design parameters. Compared with the traditional finite difference method, the sensitivity calculation time of the adjoint method is independent of the number of design variables, making it particularly suitable for the optimization problem of geometrically complex aerospace axial-flow turbomachinery blades. It significantly reduces computational resource consumption and has promising engineering application prospects.

[0022] In unsteady flow field analysis, this invention employs a harmonic balance method, transforming the unsteady flow equations into quasi-steady harmonic balance equations for solution. Compared to conventional time-domain methods, this significantly reduces computational overhead. Simultaneously, an approximate Fourier transform technique is introduced, employing a modified Gram-Schmidt orthogonalization method. By minimizing the condition number of the discrete Fourier transform and its inverse transform matrices, non-uniform time-time sampling considering multi-source disturbances is achieved. Compared to traditional oversampling techniques, this method can accurately capture the coupling effects between different disturbance sources without increasing the number of sampling points, combining high computational accuracy with low resource consumption.

[0023] In summary, the aerodynamic adjoint optimization method for aero-axial turbines of this invention fully conforms to the multi-source disturbance environment in actual aero-engine operation: it uses an adjoint method to obtain sensitivity information, a harmonic balance method to support the solution of unsteady flow fields, an approximate periodic Fourier transform method to complete non-uniform time sampling, a Hicks-Henne type function method to perform parameterized fitting of aero-axial turbine blades, and a sampling steepest descent method to update geometrically related design parameters. The greatest advantage of this invention lies in the use of a modified Gram-Schmidt orthogonalization method, which determines the non-uniform time sampling points by minimizing the condition number of the discrete Fourier transform and its inverse transform matrix. This effectively incorporates the coupling effects between multiple disturbance sources without increasing the computational burden, significantly improving the engineering applicability of the adjoint aerodynamic optimization results. Attached Figure Description

[0024] Figure 1 This is a computational grid diagram for 50% leaf height in an embodiment of the present invention.

[0025] Figure 2 This is a contour map of steady turbulent flow at 50% leaf height in an embodiment of the present invention.

[0026] Figure 3 This is the circumferential distribution and fast Fourier transform results of the 50% leaf height turbulence variables in this embodiment of the invention.

[0027] Figure 4 This is a non-uniform time sampling cloud map in an embodiment of the present invention.

[0028] Figure 5 This is a contour map of unsteady instantaneous turbulence variables at 50% of the leaf height in an embodiment of the present invention.

[0029] Figure 6 This is a comparison chart of the sensitivity of the accompanying method and the sensitivity of the finite difference method in an embodiment of the present invention.

[0030] Figure 7 This is a convergence history graph of the objective function and constraints in an embodiment of the present invention.

[0031] Figure 8 This is a comparison of the blade shape before and after optimization in an embodiment of the present invention.

[0032] Figure 9 This is a comparison of the instantaneous entropy at the outlet of an aerospace axial flow turbine before and after optimization in an embodiment of the present invention. Detailed Implementation

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

[0034] An unsteady aerodynamic adjoint optimization method for aerospace axial flow turbines considering multi-source disturbances includes the following steps; Step 1: Use software to generate fluid meshes for the aerospace axial flow turbine blades, with the mesh topology being HI monolithic structured.

[0035] Step Two: Based on the fluid computational mesh from Step One, solve the steady Navier-Stokes equations and extract the turbulence variable contour map and its circumferential distribution at 50% leaf height. The expression for the steady Navier-Stokes equations is as follows: In the formula, For flow field variables; Density; U is the velocity vector; For internal energy; For turbulence variables; F, G, and H are the convective fluxes in the axial, circumferential, and radial directions, respectively. , and , respectively, represent the viscous fluxes in the axial, circumferential, and radial directions; S is the source term; For pseudo-time; x, r and r are the axial, circumferential, and radial coordinates, respectively.

[0036] Step 3: Based on the circumferential distribution of turbulence variables at 50% blade height obtained in Step 2, perform a Fast Fourier Transform. The specific implementation process is as follows: First, the circumferential distribution of the turbulent variables at 50% blade height obtained in step two can be expressed as: In the formula, N is the number of circumferential grid points.

[0037] Secondly, regarding Performing a discrete Fourier transform yields: In the formula, It is the rotation factor; Represents the j-th Fourier coefficient; Let be the turbulence variable at the nth circumferential grid point.

[0038] Next, based on the periodicity of the Discrete Fourier Transform, the formula for the butterfly operation unit of the Fast Fourier Transform can be obtained, as shown in the following expression: In the formula, For an even number of points, use an N / 2-point DFT; For an N / 2-point DFT with an odd number of points, the expressions are as follows: Finally, in order to determine the number of harmonics required to resolve different disturbances, the following dimensionless parameters are defined: Make If this holds true, the smallest k corresponding to this condition is the number of harmonics required to resolve different disturbances.

[0039] Step 4: Based on the design speed of the aero-axial flow turbine shaft, operating conditions, and the number of blades in different blade rows, determine the unsteady flow frequencies of different disturbances and the inter-blade phase angle information corresponding to each frequency component.

[0040] The expressions for the unsteady current frequencies of the different disturbances mentioned in step four are as follows: In the formula, The frequency is unsteady. is the design rotational speed of the shaft; i is the harmonic number; k is the number of harmonics determined in step three; Operating conditions; Pi is the mathematical constant of a circle.

[0041] The expression for the inter-leaf phase angle corresponding to each frequency component is as follows: In the formula, and This indicates the number of blades in the first and second rows.

[0042] Step 5: Based on the number of harmonics of different disturbances obtained in Step 3 and the unsteady current frequency obtained in Step 4, non-uniform time point sampling is performed using the approximate periodic Fourier transform technique. The specific implementation steps of the approximate periodic Fourier transform technique described in step five for sampling at non-uniform time points are as follows: First, based on the unsteady current frequency calculated in step four, the minimum unsteady current frequency information is obtained, as shown in the following formula: Secondly, based on the minimum unsteady current frequency information, the time information is divided into M equal parts, resulting in M ​​time-sampling points. The specific calculation formula is as follows: In the formula, Let M be the sampling point at time m. In actual calculations, M is usually taken as 40000. Finally, k non-uniform sampling points are obtained using the Gram-Schmidt orthogonalization method. The specific implementation process of the Gram-Schmidt orthogonalization method is as follows: 1) Based on the time point information, construct the column vector as follows: 2) Iterative solution minimizes the condition number of the discrete Fourier transform matrix. The specific formula is as follows: In the formula, the superscript T represents the transpose of the vector; For time points The corresponding column vector; for The column vector corresponding to the given time.

[0043] 3) Obtain 2k+1 non-uniform sampling time points, as follows: In the formula, The non-uniform sampling time points are the solution to be required.

[0044] Step 6: Based on the semi-discrete form of the unsteady Navier-Stokes equations, construct the unsteady harmonic balance control equations based on the principle of harmonic balance technology, solve the unsteady harmonic balance control equations, and obtain the unsteady flow field information inside the aerospace axial flow turbine. The specific expression for the semi-discrete form of the unsteady Navier-Stokes equations is as follows: In the formula, Q represents the flow field variable, and R represents the equation space residual; Next, the flow field variables are approximated using a truncated Fourier series, expressed as: In the formula, These are time-averaged flow field variables; and denoted as Fourier coefficients of the flow field variables; k represents the number of harmonics. The interleaf phase angle; The frequency is unsteady. Then, by solving the time partial derivative of the above Fourier expression, we obtain the expression for the time spectrum source term: In the formula, For time partial derivative operators; is the Fourier coefficient vector of the flow field variables; The matrix is ​​the inverse of the discrete Fourier transform; E is the time spectrum source term operator. Finally, substituting the time spectrum source term into the unsteady flow field control equation, we obtain the unsteady harmonic balance control equation: Solving the above unsteady harmonic balance control equations yields the flow field variable Q, which is the unsteady flow field information.

[0045] Step 7: Input the unsteady flow field information described in Step 6 into the adjoint optimization solver, solve the adjoint equation, obtain the adjoint sensitivity information of the objective function with respect to the design parameters, and verify the accuracy of the adjoint sensitivity calculation using the finite difference method. The expression for the adjoint equation is as follows: In the formula, As a companion variable; I is the transpose of the time spectrum source term operator; I is the objective function, which is the mass flow rate, total pressure ratio, isentropic efficiency, or a combination thereof. This is pseudo-time.

[0046] Solving the above adjoint equation yields the adjoint variables. Substituting these adjoint variables into the sensitivity calculation formula, we can obtain the adjoint sensitivity described in step seven. The specific calculation formula is as follows: In the formula, To provide design parameters, this paper uses Hicks-Henne type functions to parameterize the blade profile. The expression for the Hicks-Henne type function is as follows: In the formula, This represents the perturbation amount of the circumferential grid coordinates; , , and The definition is as follows: In the formula, the subscripts le and te represent the leading edge and trailing edge, respectively; the subscripts hub and shroud represent the rim and hub, respectively; and M represents the radially extracted blade height number.

[0047] The expression for the objective function is: In the formula, the subscript 0 represents the initial value; This is due to entropy increase; For mass flow rate; This is the mass flow penalty function coefficient, with a value of 200; Total pressure ratio; This is the total pressure ratio penalty function coefficient, with a value of 200.

[0048] The sensitivity calculation formula for the finite difference method is as follows: In the formula, This represents the disturbance amount of the design parameters.

[0049] Step 8: Update the shape of the aerospace axial flow impeller blades based on the sensitivity information of the objective function with respect to design parameters and the steepest descent method described in Step 7; The specific expression for updating the design parameters using the steepest descent method described in step eight is as follows: In the formula, n represents the nth optimization step; This represents the disturbance step size, which is usually set to 1% of the blade chord length.

[0050] Step 9: Repeat steps 1, 6, 7 and 8 until the optimization result meets the convergence criterion, and output the optimization result.

[0051] The optimization convergence criterion described in step nine is that the L2 norm of the sensitivity is less than 1 / 2. The specific expression is as follows: To illustrate the purpose, technical solution, and advantages of the present invention in more detail, a three-stage aerospace axial flow turbine is used as an example for further detailed description. In this example, the unsteady disturbances of the middle second row of blades are respectively from the disturbances of the upstream blade row and the potential disturbances of the downstream blade row, and have multiple different disturbance sources.

[0052] Example: like Figures 1 to 9 As shown, the flow field solver and adjoint solver in this embodiment are implemented based on a self-written Fortran program, while the adjoint solver is developed based on Tapenade, an automatic differential software that converts source code. Fluid mesh generation is implemented using NUMECA's Autogrid 5, circumferential distribution extraction of turbulence variables is performed using the Data Extract function in Tecplot360, and the Fast Fourier Transform is based on the Data Transformation function in Tecplot360. The specific implementation steps of this embodiment are as follows: Step (1): The fluid computation mesh is generated using the Autogrid5 module of NUMECA. The mesh topology adopts a HI-type monolithic structured mesh and is output as a .grd file. This file contains the number of mesh nodes in the circumferential, axial, and radial directions, as well as the coordinate information of each mesh point. Figure 1 The computational mesh distribution for Blade-to-Blade fluid at 50% of the blade height is given.

[0053] Step (2): Steady flow field analysis was performed using the self-developed Fortran flow field solver TurboXD. The command ".\TurboXD -c input.txt -np 1 -flow" was executed during the calculation. Convergence was assessed by monitoring the residual curves in the res.dat file. The calculation was terminated when the residual curves became flat and the total residuals decreased by more than six orders of magnitude relative to the initial value, and flow field data files at different blade heights (e.g., 50%span.dat) were output. Extracting this file yielded a distribution cloud map of steady turbulent variables at the 50% blade height section, as shown below. Figure 2 As shown.

[0054] Step (3): Open the 50%span.data file from step (2), and use Tecplot360's Data Extract function (operation path: Data → Extract → Subzone) to extract the circumferential distribution of turbulence variables at the 50% leaf height section. Based on this, use Tecplot360's Data Transformation function to perform a Fast Fourier Transform (operation path: Data → Fourier Transformation) to obtain the Fourier coefficients corresponding to different harmonics. Finally, based on the defined dimensionless parameters... This determines the number of harmonics that need to be retained to analyze different disturbances. Figure 3 The circumferential distribution of turbulent variables at 50% leaf height and the results of fast Fourier transform are shown in the embodiments of the present invention.

[0055] Step (4): Based on the design speed, operating conditions, and number of blades in different blade rows of the aero-axial flow turbine shaft, determine the unsteady flow frequencies of different disturbances and the corresponding inter-blade phase angle information for each frequency component. The determination of the inter-blade phase angle is closely related to the number of upstream and downstream blades. Table 1 shows the unsteady flow frequencies of different blade rows and their corresponding inter-blade phase angle information.

[0056] Table 1 Step (5): Based on the number of disturbance harmonics determined in step (3) and the unsteady current frequency obtained in step (4), non-uniform time sampling is implemented using the approximate periodic Fourier transform technique. The specific operation is as follows: Modify the `time_instants` parameter in the program input file `input.txt` to `non-uniform`. After running the program, the `time_instants.txt` file is generated, which contains information about each sampling time point. In this embodiment, the distribution of non-uniform sampling time points is as follows: Figure 4 As shown.

[0057] Step (6): Set the analysis_mode parameter in the program input file input.txt to HB, and configure the number of harmonics in each row and the inter-leaf phase angle according to Table 1. Simultaneously, according to... Figure 4 The sampling time point is set based on the non-uniform sampling information. Then, the self-developed Fortran flow field solver TurboXD is run, and the input command is the same as in step (1) to solve the unsteady harmonic equilibrium equation to obtain the unsteady flow field information. Figure 5 The unsteady instantaneous turbulent variable cloud map at 50% blade height section in this embodiment is shown, demonstrating that the harmonic balance method can effectively capture the propagation characteristics of unsteady upstream and downstream flow field information between different channels.

[0058] Step (7): Using the unsteady flow field information obtained in step (6) as input, run the self-developed Fortran adjoint solver TurboADJ (execution command: .\TurboADJ -c input.txt -np 1 -adj). This solver is developed based on an automatic differentiation tool of source code conversion type. The specific process is as follows: first, use automatic differentiation software to differentiate the flow field solver, and then assemble the differentiated subroutines in reverse mode to form the adjoint solver. After solving the adjoint equation, obtain the sensitivity information of the objective function with respect to the design parameters and save it to the sensitivity.dat file. Figure 6 The comparison results of the adjoint sensitivity and the finite difference method sensitivity in this embodiment are presented. The two are in good agreement, which verifies the high accuracy of the developed adjoint solver.

[0059] Step (8): Read the sensitivity information from the sensitivity.dat file and input it into the optimizer based on the steepest descent method to update the design parameters. The optimization step size is set to 1.0e-4.

[0060] Step (9): Repeat steps (1), (6) to (8) until the optimization meets the convergence criteria. Figure 7 The convergence curves of the objective function and constraints in this embodiment are shown. After 30 optimization steps, the objective function decreases by about 2%, while the mass flow rate and total pressure ratio remain basically unchanged. Figure 8 and Figure 9 A comparison of the blade shape before and after optimization, as well as a comparison of the instantaneous entropy distribution at the outlet of the aerospace axial-flow turbine, are presented. These results further verify the reliability and effectiveness of the proposed optimization method.

Claims

1. An unsteady aerodynamic adjoint optimization method for aerospace axial-flow turbines considering multi-source disturbances, characterized in that, Includes the following steps; Step 1: Perform fluid mesh generation on the internal flow channels of the aerospace axial flow turbine, and structure the topology of the mesh generation into a single block to obtain the fluid computation mesh required for flow field analysis; Step 2: Based on the fluid computational grid obtained in Step 1, solve the steady Navier-Stokes equations and extract the turbulence variable cloud map and its circumferential distribution at 50% leaf height; Step 3: Based on the circumferential distribution of turbulent variables at 50% blade height in Step 2, perform a fast Fourier transform to determine the number of harmonics required to resolve different disturbances in the unsteady flow field analysis. Step 4: Based on the design speed, operating conditions, and number of blades in different blade rows of the aero-axial flow turbine shaft, determine the unsteady flow frequency of different disturbances and the inter-blade phase angle information corresponding to each frequency component; Step 5: Based on the number of harmonics obtained in Step 3 and the unsteady current frequency obtained in Step 4, non-uniform time point sampling is performed using the approximate periodic Fourier transform technique; Step Six: Based on the semi-discrete form of the unsteady Navier-Stokes equations, and using the number of harmonics obtained in Step Three, the unsteady flow frequency obtained in Step Four, the non-uniform time points obtained in Step Five, and the principle of harmonic balance technology, construct the unsteady harmonic balance control equations, and solve them discretely to obtain the unsteady flow field information inside the aerospace axial flow turbine. Step 7: Input the unsteady flow field information obtained in Step 6 into the adjoint optimization solver, solve the adjoint equation, obtain the adjoint sensitivity information of the objective function with respect to the design parameters, and verify the accuracy of the adjoint sensitivity calculation using the finite difference method. Step 8: Update the shape of the aerospace axial flow impeller blades based on the sensitivity information of the objective function with respect to design parameters and the steepest descent method described in Step 7; Step 9: Repeat steps 1, 6, 7 and 8 until the optimization result meets the convergence criterion; The expression for the adjoint equation described in step seven is as follows: In the formula, As a companion variable; I is the transpose of the time spectrum source term operator; I is the objective function, which is the mass flow rate, total pressure ratio, isentropic efficiency, or a combination thereof. It is pseudo-time; Solve the above adjoint equation to obtain the adjoint variables. Substitute the adjoint variables into the sensitivity calculation formula to obtain the adjoint sensitivity described in step seven. The specific calculation formula is as follows: In the formula, To provide design parameters, the blade profile is parameterized using Hicks-Henne type functions. The expression for the Hicks-Henne type function is as follows: In the formula, This represents the perturbation amount of the circumferential grid coordinates; , , and The definition is as follows: In the formula, the subscripts le and te represent the leading edge and trailing edge, respectively; the subscripts hub and shroud represent the rim and hub, respectively; M represents the radially extracted blade height number; The expression for the objective function described in step seven is: In the formula, the subscript 0 represents the initial value; This is due to entropy increase; For mass flow rate; The coefficients of the mass flow penalty function; Total pressure ratio; The total pressure ratio penalty function coefficient; The sensitivity calculation formula for the finite difference method described in step seven is as follows: In the formula, This represents the disturbance amount of the design parameters.

2. The unsteady aerodynamic adjoint optimization method for aerospace axial-flow turbines considering multi-source disturbances according to claim 1, characterized in that, The impeller blade fluid mesh generation in step one uses the Autogrid5 module of the commercial software NUMECA, and the mesh topology is HI monoblock structured.

3. The unsteady aerodynamic adjoint optimization method for aerospace axial-flow turbines considering multi-source disturbances according to claim 2, characterized in that, The expression for the steady Navier-Stokes equations in step two is as follows: In the formula, For flow field variables; Density; U is the velocity vector; For internal energy; For turbulence variables; F, G, and H are the convective fluxes in the axial, circumferential, and radial directions, respectively. , and , respectively, represent the viscous fluxes in the axial, circumferential, and radial directions; S is the source term; For pseudo-time; x, r and r are the axial, circumferential, and radial coordinates, respectively.

4. The unsteady aerodynamic adjoint optimization method for aerospace axial-flow turbofans considering multi-source disturbances according to claim 3, characterized in that, The specific implementation process of the Fast Fourier Transform described in step three is as follows: First, the circumferential distribution of the turbulent variables at 50% of the blade height obtained in step two is expressed as: In the formula, N is the number of circumferential grid points; The turbulence variable for the first grid point in the circumferential direction; For the turbulence variables at the second grid point in the circumferential direction; Let N be the turbulence variable at the Nth grid point in the circumferential direction; Secondly, regarding Performing a discrete Fourier transform yields: In the formula, It is the rotation factor; This represents the k-th Fourier coefficient; Let n be the turbulence variable at the nth circumferential grid point; Next, based on the periodicity of the Discrete Fourier Transform, the formula for the butterfly operation unit of the Fast Fourier Transform is obtained, as shown in the following expression: In the formula, the subscript j represents the order of the Fourier coefficients; N represents the number of sampling points; It is the rotation factor; Let be the j-th order Fourier coefficients; The discrete Fourier transform of N / 2 points with an even number of points; The discrete Fourier transform of N / 2 points with an odd number of points is expressed as follows: Finally, in order to determine the number of harmonics required to resolve different disturbances, the following dimensionless parameters are defined: Make If this holds true, the smallest k corresponding to this condition is the number of harmonics required to resolve different disturbances.

5. The unsteady aerodynamic adjoint optimization method for aerospace axial-flow turbines considering multi-source disturbances according to claim 4, characterized in that, The expressions for the unsteady current frequencies of the different disturbances mentioned in step four are as follows: In the formula, The frequency is unsteady. is the design rotational speed of the shaft; i is the harmonic number; Operating conditions; Pi; The expression for the inter-leaf phase angle corresponding to each frequency component described in step four is as follows: In the formula, and This indicates the number of blades in the first and second rows.

6. The unsteady aerodynamic adjoint optimization method for aerospace axial-flow turbines considering multi-source disturbances according to claim 5, characterized in that, The specific implementation steps of the approximate periodic Fourier transform technique described in step five for sampling at non-uniform time points are as follows: First, based on the unsteady current frequency calculated in step four, the minimum unsteady current frequency information is obtained, as shown in the following formula: Secondly, based on the minimum unsteady current frequency information, the time information is divided into M equal parts, resulting in M ​​time-sampling points. The specific calculation formula is as follows: In the formula, This is the sampling point at the m-th time. Finally, based on the Gram-Schmidt orthogonalization method, 2k+1 non-uniform sampling points are obtained. The specific implementation process of the Gram-Schmidt orthogonalization method is as follows: 1) Based on the time point information, construct the column vector as follows: 2) Iterative solution minimizes the condition number of the discrete Fourier transform matrix. The specific formula is as follows: In the formula, the superscript T denotes the transpose of the vector; For time points The corresponding column vector; for The column vector corresponding to the given time; 3) Obtain 2k+1 non-uniform sampling time points, as follows: In the formula, The non-uniform sampling time points are the solution to be required.

7. The unsteady aerodynamic adjoint optimization method for aerospace axial-flow turbofans considering multi-source disturbances according to claim 6, characterized in that, The specific expression for the semi-discrete form of the unsteady Navier-Stokes equations described in step six is ​​as follows: In the formula, Q represents the flow field variable, and R represents the equation space residual; The flow field variables are approximated using truncated Fourier series, expressed as follows: In the formula, These are time-averaged flow field variables; and These are the Fourier coefficients of the flow field variables; The interleaf phase angle; The frequency is unsteady. For time; Pi; Solving for the time partial derivative of the above Fourier expression yields the expression for the time spectrum source term, as follows: In the formula, For time partial derivative operators; is the Fourier coefficient vector of the flow field variables; The matrix is ​​the inverse of the discrete Fourier transform; E is the time spectrum source term operator. Finally, substituting the time spectrum source term into the unsteady flow field control equation, we obtain the aforementioned unsteady harmonic balance control equation: Solving the above unsteady harmonic equilibrium control equations yields the flow field variable Q, thus providing the unsteady flow field information.

8. The unsteady aerodynamic adjoint optimization method for aerospace axial-flow turbines considering multi-source disturbances according to claim 7, characterized in that, The specific expression for updating the design parameters using the steepest descent method described in step eight is as follows: In the formula, n represents the nth optimization step; This indicates the perturbation step size.

9. The unsteady aerodynamic adjoint optimization method for aerospace axial-flow turbines considering multi-source disturbances according to claim 8, characterized in that, Step nine optimizes the convergence criterion by setting the sensitivity L2 norm to less than 1. The specific expression is as follows: 。

Images