A multi-objective optimization method for aerodynamic profile of full-scale subsonic wind tunnel loop based on steady CFD

Abstract

The application provides a full-size subsonic wind tunnel loop aerodynamic profile multi-objective optimization method based on steady CFD, and belongs to the technical field of wind tunnels.The application reduces the initial geometric parameters to a sensitive parameter subspace through principal component analysis, generates initial sample points by using Latin hypercube sampling, extracts response values by executing steady Reynolds average Navier-Stokes equation solving, constructs a Kriging surrogate model, performs multi-start parallel search based on an expected improvement function, extracts a separation zone length and a backflow zone height to input an unsteady correction factor model to calculate a correction coefficient for correcting a diffusion section pressure recovery coefficient, searches a Pareto front solution set by using a non-dominated sorting genetic algorithm, selects an optimal parameter combination through a comprehensive performance index, and reconstructs the wind tunnel loop aerodynamic profile through a hierarchical parameterization method and performs verification calculation, so that the problem that the steady calculation framework cannot accurately predict the separation flow and the optimization result deviates from the true performance is solved.

Description

Technical Field

[0001] This invention belongs to the field of wind tunnel technology, and more specifically, relates to a multi-objective optimization method for the aerodynamic profile of a full-size subsonic wind tunnel loop based on steady CFD. Background Technology

[0002] Optimizing the aerodynamic profile of a full-scale subsonic wind tunnel loop requires coordinating multiple performance indicators, such as pressure recovery in the diffuser section, velocity uniformity in the contraction section, and total pressure loss in the corner region. Traditional methods use computational fluid dynamics (CFD) to evaluate the aerodynamic performance of different geometries and obtain the optimal profile design through iterative search. In current wind tunnel loop optimization practices, due to the extremely high computational cost of unsteady Reynolds-averaged Navier-Stokes equations or separated vortex simulations, engineering optimization mainly relies on steady Reynolds-averaged Navier-Stokes equation solvers for performance evaluation. However, steady computational frameworks, based on turbulence model assumptions, struggle to accurately capture flow separation and backflow phenomena within the diffuser section. This leads to systematic biases in predicting the separation zone scale and pressure recovery efficiency, causing optimization algorithms to converge to a geometric parameter combination that is not optimal for the actual flow field, based on distorted objective function searches. Consequently, the optimization results deteriorate significantly under actual unsteady flow conditions. In other words, existing technologies suffer from the technical problem that steady computational frameworks cannot accurately predict separated flows, resulting in optimization results deviating from actual performance. Summary of the Invention

[0003] In view of this, the present invention provides a multi-objective optimization method for the aerodynamic profile of a full-size subsonic wind tunnel loop based on steady CFD, which can solve the technical problem in the prior art where the steady calculation framework cannot accurately predict the separation flow in the multi-objective optimization process of the aerodynamic profile of a full-size subsonic wind tunnel loop, resulting in the optimization results deviating from the actual performance.

[0004] This invention is implemented as follows: It provides a multi-objective optimization method for the aerodynamic profile of a full-size subsonic wind tunnel loop based on steady-state CFD, comprising: acquiring initial geometric parameters and setting boundary condition parameters; performing principal component analysis on the initial geometric parameters to form a dimension-reduced parameter vector; generating an initial sample point set in the dimension-reduced parameter vector space using Latin hypercube sampling and extracting response values ​​by solving the steady-state Reynolds-averaged Navier-Stokes equations; constructing a Kriging surrogate model; calculating the expected improvement function value based on the Kriging surrogate model and performing a multi-starting-point parallel search to extract the separation zone length and recirculation zone height; inputting the separation zone length and recirculation zone height into an unsteady correction factor model to calculate the correction coefficient and obtain the corrected diffusion section pressure recovery coefficient; constructing a multi-objective optimization function and using a non-dominated sorting genetic algorithm to search for the Pareto front solution set; calculating the comprehensive performance index, selecting the optimal parameter combination, and reconstructing the aerodynamic profile of the wind tunnel loop; performing verification calculations and outputting the optimization results.

[0005] The initial geometric parameters are specifically the length of the diffuser section, the outlet area of ​​the diffuser section, the corner radius, and the coordinates of the control points of the contraction section profile. The boundary condition parameters are specifically the total inlet pressure, the total inlet temperature, and the static outlet pressure.

[0006] Specifically, the principal component analysis steps involve constructing the covariance matrix of the initial geometric parameters and solving for eigenvalues ​​and eigenvectors to achieve dimensionality reduction, identifying principal component directions where the ratio of eigenvalues ​​to the sum of eigenvalues ​​is greater than 0.85, and projecting the initial geometric parameters onto the sensitive parameter subspace.

[0007] The sensitive parameter subspace refers to the low-dimensional parameter space that has a significant impact on the optimization objective, which is screened out through principal component analysis, reducing the search dimension of the optimization problem from the original 20 to 30 dimensions to 3 to 5 dimensions.

[0008] The Latin hypercube sampling mentioned above refers to a stratified sampling technique that divides the range of values ​​for each parameter into several intervals and extracts exactly one sample from each interval.

[0009] Specifically, the step of extracting response values ​​by solving the steady Reynolds-averaged Navier-Stokes equation involves constructing a three-dimensional computational domain grid for each sample point in the initial sample point set and solving the steady Reynolds-averaged Navier-Stokes equation to extract the uniformity of the exit velocity in the contraction section, the pressure recovery coefficient in the diffusion section, and the total pressure loss coefficient in the corner region as response values.

[0010] The Kriging agent model refers to an alternative model based on spatial interpolation theory, which assumes that the response value is composed of a deterministic trend function and a stochastic process, and determines the hyperparameters of the correlation function through maximum likelihood estimation.

[0011] Specifically, the step of constructing the Kriging surrogate model involves using the dimensionality-reduced parameter vector and response values ​​of the initial sample point set to construct the Kriging surrogate model, calculating the root mean square error between the predicted value of the Kriging surrogate model and the calculated value of the computational fluid dynamics, and proceeding to the next step when the root mean square error is less than 5%; otherwise, additional sample points are added and the Kriging surrogate model is updated.

[0012] Specifically, the expected improvement function value is determined by comprehensively considering the gap between the Kriging surrogate model prediction value and the current optimal objective function value, as well as the prediction uncertainty, and selecting the top 20 candidate points with the largest expected improvement function values ​​as the initial positions for multi-starting point parallel search.

[0013] Specifically, the step of performing a multi-starting-point parallel search to extract the separation zone length and the recirculation zone height involves solving the steady Reynolds-averaged Navier-Stokes equation at each initial position and extracting the separation zone length and the recirculation zone height. The separation zone length and the recirculation zone height are used as flow separation characteristic quantities under the steady calculation framework.

[0014] Specifically, the unsteady correction factor model is established by calculating a small number of unsteady Reynolds-averaged Navier-Stokes equations or by vortex simulation to establish an empirical correction relationship. The separation zone length and recirculation zone height are used as input variables. The functional relationship between the correction coefficient and the separation zone length and recirculation zone height is fitted by a neural network or polynomial regression. The separation Reynolds number and the reverse pressure gradient parameters are calculated, and the correction coefficient is calculated based on the separation Reynolds number and the reverse pressure gradient parameters.

[0015] The separation Reynolds number refers to a dimensionless parameter characterizing the flow separation tendency, defined as the product of the boundary layer momentum thickness and the incoming flow velocity at the separation point divided by the kinematic viscosity. The adverse pressure gradient parameter is specifically calculated by dividing the friction static pressure difference by the product of the dynamic pressure and the characteristic length.

[0016] Specifically, the step of obtaining the corrected diffusion section pressure recovery coefficient involves multiplying the correction coefficient by the diffusion section pressure recovery coefficient obtained by solving the steady Reynolds-averaged Navier-Stokes equations to obtain the corrected diffusion section pressure recovery coefficient. The step of constructing the multi-objective optimization function involves obtaining the corrected diffusion section pressure recovery coefficient, the uniformity of the contraction section outlet velocity, and the total pressure loss coefficient in the corner region, and then normalizing them by dividing them by the corresponding reference values ​​to construct the multi-objective optimization function.

[0017] Specifically, the non-dominated sorting genetic algorithm sorts the population individuals in a hierarchical manner according to their dominance relationship. Non-dominated individuals have higher fitness, and the diversity of the solution set distribution is maintained by crowding distance. The Pareto front solution set is searched in the dimensionality-reduced parameter vector space.

[0018] Specifically, the comprehensive performance index is obtained by multiplying the normalized uniformity of the contraction section outlet velocity by a weighting coefficient of 0.35, the normalized corrected pressure recovery coefficient of the diffusion section by a weighting coefficient of 0.45, and the reciprocal of the normalized total pressure loss coefficient of the corner region by a weighting coefficient of 0.20, and then summing the results. The parameter combination with the highest comprehensive performance index is selected.

[0019] The step of reconstructing the aerodynamic profile of the wind tunnel loop specifically involves inversely transforming the parameter combination from the sensitive parameter subspace to the complete geometric parameter space, and reconstructing the aerodynamic profile of the wind tunnel loop using a hierarchical parameterization method. This hierarchical parameterization method specifically employs a multi-scale geometric expression strategy. The large-scale profile uses 6 to 8 B-spline curve control points to define the overall orientation of the diffuser and contraction sections, while small-scale local features are superimposed with a free-deformable mesh. The step of performing verification calculations specifically involves performing steady Reynolds-averaged Navier-Stokes equations verification calculations on the reconstructed aerodynamic profile of the wind tunnel loop. When the uniformity of the exit velocity of the contraction section in the verification calculation is greater than 95% and the pressure recovery coefficient of the corrected diffuser section in the verification calculation is greater than 0.75, the reconstructed aerodynamic profile of the wind tunnel loop is output as the optimization result.

[0020] This invention achieves parameter space dimensionality reduction through principal component analysis and reduces the number of computational fluid dynamics calls by combining a Kriging surrogate model. Within the steady-state calculation framework, the length of the separation zone and the height of the recirculation zone are extracted as flow separation characteristics. An unsteady correction factor model is constructed to establish the mapping relationship between these characteristics and the actual pressure recovery coefficient. The correction factor is applied to the steady-state calculation results to compensate for prediction biases in the separation flow. The unsteady correction factor model establishes empirical correction relationships through a small number of high-precision unsteady calculations or separation vortex simulations. This allows the optimization process to achieve performance evaluations close to the unsteady real flow field while maintaining the efficiency of steady-state calculations, avoiding convergence to a pseudo-optimal solution due to objective function distortion. This ensures that the final output wind tunnel loop aerodynamic profile still meets the engineering threshold requirements for the uniformity of the exit velocity in the contraction section and the pressure recovery coefficient in the diffusion section under real separation flow conditions. In summary, this invention solves the technical problem mentioned in the background art where the steady-state calculation framework cannot accurately predict separation flow, leading to optimization results deviating from actual performance. Attached Figure Description

[0021] Figure 1 This is a flowchart of the method of the present invention.

[0022] Figure 2 This is a diagram showing the distribution of the length of the diffusion section separation zone and the height of the recirculation zone in the embodiment.

[0023] Figure 3 This is a comparison chart of the steady-state calculation and the corrected pressure recovery coefficient of the diffusion section in the embodiment. Detailed Implementation

[0024] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below.

[0025] The invention provides a multi-objective optimization method for the aerodynamic profile of a full-scale subsonic wind tunnel loop based on steady CFD, including:

[0026] S01. Collect the initial geometric parameters of the full-size subsonic wind tunnel loop. The initial geometric parameters include the length of the diffuser section, the outlet area of ​​the diffuser section, the corner radius, and the coordinates of the control points of the contraction section profile. Set the boundary condition parameters, including the inlet total pressure, the inlet total temperature, and the outlet static pressure.

[0027] S02. Perform principal component analysis on the initial geometric parameters, calculate the eigenvalues ​​of the covariance matrix of the initial geometric parameters, identify the principal component directions whose ratio of the eigenvalue to the sum of the eigenvalues ​​is greater than 0.85, and project the initial geometric parameters to the sensitive parameter subspace to form a dimension-reduced parameter vector.

[0028] S03. In the reduced-dimensional parameter vector space, an initial sample point set is generated by Latin hypercube sampling. A three-dimensional computational domain grid is constructed for each sample point in the initial sample point set, and the steady Reynolds-averaged Navier-Stokes equation is solved. The uniformity of the outlet velocity of the contraction section, the pressure recovery coefficient of the diffusion section, and the total pressure loss coefficient of the corner region are extracted as response values.

[0029] S04. Construct a Kriging surrogate model using the dimensionality-reduced parameter vector and response value of the initial sample point set. Calculate the root mean square error between the predicted value of the Kriging surrogate model and the calculated value of the computational fluid dynamics. If the root mean square error is less than 5%, proceed to step S05. Otherwise, add supplementary sample points and update the Kriging surrogate model, then repeat step S04.

[0030] S05. Calculate the expected improvement function value according to the Kriging surrogate model, select the top 20 candidate points with the largest expected improvement function value as the initial positions of the multi-starting point parallel search, and perform steady Reynolds-averaged Navier-Stokes equations to solve each initial position and extract the separation zone length and recirculation zone height.

[0031] S06. Input the separation zone length and recirculation zone height into the unsteady correction factor model, calculate the separation Reynolds number and reverse pressure gradient parameters through the unsteady correction factor model, calculate the correction coefficient based on the separation Reynolds number and reverse pressure gradient parameters, and multiply the correction coefficient by the diffusion section pressure recovery coefficient obtained by solving the steady Reynolds-averaged Navier-Stokes equation in step S05 to obtain the corrected diffusion section pressure recovery coefficient.

[0032] S07. Obtain the corrected diffusion section pressure recovery coefficient, the contraction section outlet velocity uniformity extracted in step S05, and the corner region total pressure loss coefficient extracted in step S05. After normalizing the corrected diffusion section pressure recovery coefficient, the contraction section outlet velocity uniformity, and the corner region total pressure loss coefficient by the corresponding reference values, construct a multi-objective optimization function. Use a non-dominated sorting genetic algorithm to search for the Pareto front solution set in the reduced-dimensional parameter vector space.

[0033] S08. Calculate the comprehensive performance index of each solution in the Pareto front solution set. The comprehensive performance index is obtained by multiplying the normalized contraction section exit velocity uniformity by a weighting coefficient of 0.35, the normalized corrected diffusion section pressure recovery coefficient by a weighting coefficient of 0.45, and the reciprocal of the normalized corner region total pressure loss coefficient by a weighting coefficient of 0.20, and then summing them up. Select the parameter combination with the highest comprehensive performance index, and inversely transform the parameter combination from the sensitive parameter subspace to the complete geometric parameter space. Reconstruct the wind tunnel loop aerodynamic profile using a hierarchical parameterization method.

[0034] S09. Perform steady Reynolds-averaged Navier-Stokes equation verification calculations on the reconstructed wind tunnel loop aerodynamic profile, obtain the uniformity of the exit velocity of the contraction section and the pressure recovery coefficient of the corrected diffusion section from the verification calculation. When the uniformity of the exit velocity of the contraction section from the verification calculation is greater than the set threshold for the uniformity of the exit velocity of the contraction section and the pressure recovery coefficient of the corrected diffusion section from the verification calculation is greater than the set threshold for the pressure recovery coefficient of the diffusion section, output the reconstructed wind tunnel loop aerodynamic profile as the optimization result; otherwise, adjust the separation Reynolds number correction weight and the inverse pressure gradient parameter correction weight in the unsteady correction factor model and return to step S06.

[0035] The threshold value for the uniformity of the exit velocity in the contraction section is set to 95%, and the threshold value for the pressure recovery coefficient in the diffusion section is set to 0.75. The weighting coefficients of 0.35, 0.45, and 0.20 are calculated using the analytic hierarchy process (AHP). The AHP constructs a judgment matrix based on the wind tunnel loop performance requirements. In the judgment matrix, the importance ratio of the uniformity of the exit velocity in the contraction section to the pressure recovery coefficient in the diffusion section is 1:1.3, the importance ratio of the uniformity of the exit velocity in the contraction section to the total pressure loss coefficient in the corner region is 1.75:1, and the importance ratio of the pressure recovery coefficient in the diffusion section to the total pressure loss coefficient in the corner region is 2.25:1. The weighting coefficients are obtained by normalizing the eigenvectors of the judgment matrix.

[0036] Principal component analysis (PCA) is a dimensionality reduction method that constructs a covariance matrix of initial geometric parameters and solves for eigenvalues ​​and eigenvectors. Each element of the covariance matrix represents the degree of linear correlation between two geometric parameters, the eigenvalues ​​reflect the variance contribution rate of the corresponding principal component direction, and the eigenvectors define the new coordinate axis directions in the parameter space. When the ratio of the eigenvalue of a principal component to the sum of all eigenvalues ​​exceeds 0.85, the principal component is considered to contain the main information of the original parameters, and other principal components with smaller eigenvalues ​​are discarded. This compresses the high-dimensional space containing dozens of parameters into a low-dimensional sensitive parameter subspace composed of 3 to 5 principal components.

[0037] The sensitive parameter subspace is a low-dimensional parameter space selected through principal component analysis that significantly impacts the optimization objective. In the original geometric parameter space, the influence of parameters such as the diffuser length and corner radius on flow field characteristics varies significantly. Small changes in some parameters can lead to drastic fluctuations in the pressure recovery coefficient of the diffuser or the uniformity of the outlet velocity of the contraction section, while adjustments to other parameters have almost no effect on the objective function. The sensitive parameter subspace retains only the parameter directions corresponding to the former, reducing the search dimension of the optimization problem from the original 20-30 dimensions to 3-5 dimensions, significantly reducing the number of computational fluid dynamics calculations while maintaining optimization accuracy.

[0038] Latin hypercube sampling is a stratified sampling technique that divides the range of each parameter into several equal intervals, ensuring that exactly one sample is drawn from each interval, resulting in a uniform distribution of sample points within the parameter space. Compared to random sampling, Latin hypercube sampling covers a wider parameter space with fewer samples, improving the global approximation accuracy of the surrogate model.

[0039] The Kriging surrogate model is an alternative model based on spatial interpolation theory, used to predict response values ​​at unknown design points under finite sample conditions. The Kriging surrogate model assumes that the response value is composed of a deterministic trend function and a stochastic process. It determines the hyperparameters of the correlation function through maximum likelihood estimation, ensuring that the model's predicted values ​​accurately match the true values ​​at known sample points and provide the optimal unbiased estimate at unknown points. The Kriging surrogate model not only outputs the predicted values ​​but also provides the prediction variance as a measure of uncertainty, providing necessary information for calculating the expected improvement criteria.

[0040] The expected improvement function value is an indicator that measures the potential improvement of a candidate design point relative to the current optimal solution. This expected improvement function value comprehensively considers the difference between the Kriging surrogate model prediction and the current optimal objective function value, as well as prediction uncertainty. When the predicted value of a design point is significantly better than the current optimal value and the prediction variance is large, the expected improvement function value is high, indicating that the design point is worthwhile for expensive computational fluid dynamics verification calculations. The expected improvement criterion balances the need to develop known optimal regions with the need to explore unknown regions, avoiding premature convergence of the optimization algorithm to a local optimum.

[0041] The multi-starting-point parallel search strategy involves simultaneously launching the optimization algorithm from multiple different initial positions. By exploring different regions of the parameter space in parallel, the risk of getting trapped in local optima is reduced. The top 20 candidate points with the largest expected improvement function values ​​are selected as initial positions to ensure that the search covers the potential region of the global optimum. Each starting point independently performs the solution of the steady Reynolds-averaged Navier-Stokes equations. Finally, the convergence results of each starting point are compared and the optimal solution is selected.

[0042] The unsteady correction factor model is an empirical correction relationship established through calculations using a small number of unsteady Reynolds-averaged Navier-Stokes equations or separated eddy simulations. This unsteady correction factor model uses the separation zone length and recirculation zone height, which cannot be accurately predicted in steady-state calculations, as input variables. It then uses neural networks or polynomial regression to fit the functional relationship between the correction coefficient and these dimensions. Multiplying the correction coefficient by the pressure recovery coefficient of the diffuser section from the steady-state calculations compensates for the systematic bias caused by the insufficient prediction of separated flows by the steady-state Reynolds-averaged Navier-Stokes equations. This allows the optimization process to achieve performance evaluation results close to the actual unsteady flow field within the steady-state calculation framework.

[0043] The separation Reynolds number is a dimensionless parameter characterizing the tendency of flow separation. It is defined as the product of the boundary layer momentum thickness and the incoming flow velocity at the separation point, divided by the kinematic viscosity. A higher separation Reynolds number indicates a stronger resistance of the boundary layer to the adverse pressure gradient, making flow separation more difficult. Within the diffuser section, as the cross-sectional area increases, the static pressure rises, leading to an enhanced adverse pressure gradient. When the separation Reynolds number falls below a critical value, boundary layer separation occurs, forming a large-scale recirculation region that reduces the pressure recovery efficiency of the diffuser section.

[0044] The adverse pressure gradient parameter is a dimensionless quantity describing the rate of increase in static pressure along the flow direction. It is calculated by dividing the static pressure difference along the flow path by the product of the dynamic pressure and the characteristic length. In the diffuser section of the wind tunnel, the expansion of the cross-sectional area leads to a decrease in flow velocity and an increase in static pressure, forming an unfavorable adverse pressure gradient. The adverse pressure gradient parameter quantifies the intensity of this reverse pressure increase. When the adverse pressure gradient parameter exceeds 0.02, the low-momentum fluid in the boundary layer experiences backflow under pressure, triggering flow separation.

[0045] Among them, the non-dominated sorting genetic algorithm is an evolutionary algorithm for handling multi-objective optimization problems. It searches for Pareto optimal solutions by simulating natural selection and genetic mechanisms. The non-dominated sorting genetic algorithm hierarchically sorts individuals in the population according to their dominance relationship. Non-dominated individuals have higher fitness, and the algorithm maintains the diversity of the solution distribution through crowding distance. After iterative evolution through selection, crossover, and mutation genetic operations, the population gradually approaches the Pareto front, outputting a set of mutually non-dominated optimal solutions. This provides decision-makers with a trade-off between the pressure recovery coefficient of the diffusion section, the uniformity of the exit velocity of the contraction section, and the total pressure loss coefficient in the corner region.

[0046] The Pareto front solution set is the set of all Pareto optimal solutions in a multi-objective optimization problem. A solution is called Pareto optimal if and only if there is no other solution that is not inferior to the solution in all objectives and is superior to the solution in at least one objective. In wind tunnel loop optimization, increasing the pressure recovery coefficient in the diffuser section often comes at the cost of increasing the total pressure loss in the corner region, while improving the uniformity of the exit velocity in the contraction section limits the diffuser angle, leading to an increase in loop length. The Pareto front solution set demonstrates the irreconcilable conflict between these objectives, with each solution representing a different performance balance scheme.

[0047] The Analytic Hierarchy Process (AHP) is a multi-criteria decision-making method that combines quantitative and qualitative approaches. It constructs a judgment matrix to express the relative importance of different objectives. The elements of the judgment matrix represent the importance ratios of pairwise comparisons between objectives. By calculating the eigenvector corresponding to the largest eigenvalue of the judgment matrix and normalizing it, the weight coefficients of each objective are obtained. These weight coefficients reflect the differences in the degree of importance placed on the uniformity of the exit velocity in the contraction section, the pressure recovery coefficient in the diffusion section, and the total pressure loss coefficient in the corner region during engineering applications.

[0048] The hierarchical parameterization method employs a multi-scale geometric representation strategy to reconstruct the wind tunnel loop profile. The large-scale profile defines the overall orientation of the diffuser and contraction sections using 6 to 8 B-spline curve control points, ensuring surface continuity and smoothness. Small-scale local features are superimposed with a free-deformable mesh, applying local curvature adjustments at corner transition zones and the beginning of the diffuser section to capture detailed geometry sensitive to flow separation. Shape syntax constraints ensure the reconstructed profile meets curvature continuity requirements and is fabricable; control point displacement is limited to ±15% of the initial position to prevent unreasonable abrupt changes in elevation.

[0049] Among them, the separation Reynolds number correction weight and the adverse pressure gradient parameter correction weight are coefficients used in the unsteady correction factor model to adjust the contribution ratio of the separation Reynolds number and the adverse pressure gradient parameter to the correction coefficient. When the verification calculation results do not meet the set threshold for the uniformity of the exit velocity in the contraction section or the set threshold for the pressure recovery coefficient in the diffusion section, the correction strength for boundary layer separation is enhanced by increasing the separation Reynolds number correction weight, or the correction strength for pressure recovery is enhanced by increasing the adverse pressure gradient parameter correction weight, so that the corrected pressure recovery coefficient in the diffusion section is closer to the actual unsteady flow results.

[0050] The 95% threshold for the uniformity of the exit velocity in the contraction section is derived from the engineering requirements for flow field quality in the full-scale subsonic wind tunnel test section. When the uniformity of the exit velocity in the contraction section is lower than 95%, the uneven velocity distribution within the test section will cause the measurement error to exceed the allowable range. The 0.75 threshold for the pressure recovery coefficient in the diffuser section is derived from the wind tunnel loop energy loss control standard. When the pressure recovery coefficient in the diffuser section is lower than 0.75, the wind tunnel drive power consumption exceeds the economic operating limit.

[0051] The present invention also provides a multi-objective optimization system for the aerodynamic profile of a full-size subsonic wind tunnel loop based on steady CFD, implemented by a computer. The computer is equipped with a storage medium, which stores program instructions. When the program instructions are run in the computer, they execute the above-described method.

[0052] The specific implementation methods of the above steps are described in detail below.

[0053] The specific implementation of step S01 is to obtain the original geometric data of the wind tunnel loop through measuring instruments or computer-aided design software, including the axial length of the diffuser section, the cross-sectional area of ​​the outlet, the radius of the arc at the corner, and the three-dimensional coordinates of the control points of the Bézier curve of the contraction section. At the same time, according to the wind tunnel design specifications, the inlet total pressure is set to 1.05 times the standard atmospheric pressure, the inlet total temperature is 288K, and the outlet static pressure is the standard atmospheric pressure. These parameters constitute the complete input conditions for subsequent computational fluid dynamics simulation, ensuring that the optimization process is based on real and feasible engineering constraints.

[0054] The specific implementation of step S02 is as follows: First, the collected initial geometric parameters are constructed into a matrix. The covariance matrix between each column of the matrix is ​​calculated to quantify the correlation between parameters. The covariance matrix is ​​decomposed into eigenvalues ​​to obtain eigenvalues ​​and corresponding eigenvectors. After sorting the eigenvalues ​​from largest to smallest, the ratio of each eigenvalue to the sum of all eigenvalues ​​is calculated. The top principal components with a cumulative contribution rate of more than 85% are selected as sensitive directions. The original high-dimensional parameter vector is multiplied by the eigenvector matrix corresponding to these principal components to complete the projection transformation, resulting in a dimensionality-reduced parameter vector with dimensions of 3 to 5. This process utilizes the orthogonal transformation principle in linear algebra to compress the parameter space from 20 to 30 dimensions to a low-dimensional subspace, significantly reducing the dimensionality burden of subsequent optimization calculations.

[0055] The specific implementation of step S03 involves applying Latin hypercube sampling technology to generate initial sample points in the dimensionality-reduced parameter space. This technology divides the value range of each parameter into several equally spaced sub-intervals and randomly selects a sample value within each sub-interval. Through permutation and combination, it ensures that the sample points are uniformly distributed in the multidimensional space. For the geometric configuration corresponding to each sample point, computational fluid dynamics preprocessing software is used to generate a three-dimensional computational domain containing boundary layer meshes and volume meshes. The total number of meshes is controlled within a certain range. to The magnitude is adjusted to balance computational accuracy and efficiency. Then, the steady Reynolds-averaged Navier-Stokes equation solver is called to calculate the flow field. From the converged flow field results, the ratio of the standard deviation of the velocity distribution at the exit section of the contraction section to the average velocity is extracted as the velocity uniformity, the ratio of the static pressure difference between the outlet and inlet of the diffuser section to the dynamic pressure is extracted as the pressure recovery coefficient, and the ratio of the total pressure loss in the corner region to the inlet dynamic pressure is extracted as the total pressure loss coefficient. These response values ​​directly reflect the aerodynamic performance of the wind tunnel loop.

[0056] The specific implementation of step S04 involves using the dimensionality-reduced parameter vector of existing sample points as input and the corresponding three response values ​​as output. A spatial substitution model is constructed using Kriging interpolation theory. This model assumes that the response value is composed of a global trend term and a local correlation random term. The length scale parameter and smoothness parameter of the correlation function are determined by the maximum likelihood estimation method, so that the predicted value of the model at the known sample points is completely consistent with the actual calculated value. The model is cross-validated to calculate the root mean square error between the predicted value and the computational fluid dynamics calculated value. When the error exceeds 5%, it indicates that the sample point coverage is insufficient. At this time, it is necessary to add supplementary sample points in the region with the largest model prediction variance and retrain the model. This process is iterated until the model accuracy meets the requirements. The purpose of this step is to establish a low-cost approximate model to replace expensive flow field simulation.

[0057] The specific implementation of step S05 is to calculate the expected improvement function based on the trained Kriging surrogate model. This function comprehensively considers the probability that the predicted target value at a certain point is better than the current optimal value and the magnitude of the prediction uncertainty at that point. By evaluating the expected improvement function value in the entire parameter space, the candidate points with the most potential for improvement are identified. The top 20 points with the highest expected improvement function values ​​are selected as the initial positions for multi-starting point optimization. The steady Reynolds-averaged Navier-Stokes equations are solved for the geometric configuration corresponding to each initial position. In the convergent flow field, the flow direction length of the separation zone is determined by identifying the isosurface where the axial velocity is zero, and the normal height of the recirculation zone is determined by identifying the location of the maximum reverse velocity in the recirculation zone. These two parameters characterize the severity of flow separation. The multi-starting point parallel search strategy can avoid the optimization algorithm from getting trapped in local optima.

[0058] The specific implementation of step S06 involves inputting the extracted separation zone length and recirculation zone height into a pre-established unsteady correction factor model. This model is obtained through neural network training based on a small number of high-precision unsteady Reynolds-averaged Navier-Stokes equations or separation vortex simulation results. Internally, the model first calculates the separation Reynolds number as the boundary layer momentum thickness at the separation point multiplied by the incoming velocity and then divided by the fluid kinematic viscosity. Simultaneously, it calculates the adverse pressure gradient parameter as the friction static pressure difference divided by the product of dynamic pressure and characteristic length. After multiplying these two dimensionless parameters by their respective correction weights, the correction coefficient is obtained through nonlinear function mapping. This correction coefficient reflects the systematic deviation of steady calculations in predicting separated flows. Multiplying the correction coefficient by the diffusion section pressure recovery coefficient obtained from steady calculations in step S05 yields a corrected pressure recovery coefficient that closely approximates the actual unsteady flow. The core function of this step is to compensate for the insufficient prediction of complex separated flows by the steady model while maintaining the efficiency of steady calculations.

[0059] The specific implementation of step S07 involves obtaining three optimization objectives: the corrected diffusion section pressure recovery coefficient, the contraction section outlet velocity uniformity extracted in step S05, and the total pressure loss coefficient in the corner region. These objectives are then normalized by dividing them by reference values ​​of 0.70, 0.92, and 0.05, respectively, to ensure that each objective is of the same magnitude. A weighted comprehensive objective function is constructed as follows: normalized velocity uniformity multiplied by 0.35, plus normalized pressure recovery coefficient multiplied by 0.45, minus normalized total pressure loss coefficient multiplied by 0.20. The weighting coefficients are determined using the analytic hierarchy process (AHP), based on the importance ratio of velocity uniformity to pressure recovery coefficient (1:1.3). A judgment matrix is ​​constructed using a total pressure loss coefficient importance ratio of 1.75:1 and a pressure recovery coefficient importance ratio of 2.25:1. The eigenvector corresponding to the largest eigenvalue of this matrix is ​​calculated and normalized to obtain the weights. Subsequently, a non-dominated sorting genetic algorithm is applied to perform multi-objective optimization in the dimensionality-reduced parameter space. This algorithm maintains a population of 100 to 200 individuals, sorts individuals in a hierarchical manner through non-dominated relations, and calculates crowding distance to maintain the diversity of solutions. After 50 to 100 generations of genetic operations such as selection, crossover, and mutation, it converges to the Pareto front solution set, which contains a series of optimal solutions that achieve different balances among the three objectives.

[0060] The specific implementation of step S08 involves calculating a comprehensive performance index for each solution in the Pareto front solution set. This index is obtained by multiplying the normalized velocity uniformity corresponding to the solution by 0.35, the normalized pressure recovery coefficient after correction by 0.45, and the reciprocal of the normalized total pressure loss coefficient by 0.20, and then summing the results. The parameter combination with the highest comprehensive performance index value is selected as the final optimization result. This parameter combination is then mapped from the sensitive parameter subspace back to the complete geometric parameter space through principal component inverse transformation to recover the specific values ​​of all geometric parameters. A hierarchical parameterization method is used to reconstruct the aerodynamic profile of the wind tunnel loop. This method first uses 6 to 8 Bézier spline curve control points to define the overall profile of the diffusion and contraction sections to ensure the global smoothness of the surface. Then, a free-deformable mesh is superimposed in the corner transition area to adjust the local curvature in order to capture the detailed features sensitive to flow separation. During the reconstruction process, geometric constraints are applied to ensure that the curvature of adjacent surfaces is continuous and the displacement of the control points does not exceed ±15% of the initial position to ensure manufacturing feasibility. The purpose of this step is to transform the abstract parameters obtained from the optimization into specific geometric profiles that can be directly used for engineering implementation.

[0061] The specific implementation of step S09 involves recalculating the reconstructed wind tunnel loop aerodynamic profile using the steady Reynolds-averaged Navier-Stokes equations to verify it. The velocity uniformity at the exit of the contraction section and the pressure recovery coefficient of the corrected diffusion section are extracted from the verification calculation. It is determined whether the velocity uniformity is greater than 95% and the pressure recovery coefficient is greater than 0.75. If both conditions are met, the optimization result meets the wind tunnel performance requirements, and the reconstructed profile is output as the final solution. If either condition is not met, the weights for the separation Reynolds number correction and the reverse pressure gradient parameter correction in the unsteady correction factor model need to be adjusted. When the pressure recovery coefficient is low, the weight of the reverse pressure gradient parameter correction is increased to strengthen the correction of pressure recovery. When the velocity uniformity is insufficient, the weight of the separation Reynolds number correction is increased to strengthen the correction of boundary layer separation. After adjustment, the process returns to step S06 to recalculate the corrected performance indicators. This iterative process ensures that the performance evaluation obtained under the steady calculation framework of the final optimization result matches the actual unsteady flow characteristics. The velocity uniformity threshold of 95% originates from the flow field quality requirements of the test section, and the pressure recovery coefficient threshold of 0.75 originates from the wind tunnel energy consumption economic constraints.

[0062] It should be noted that the first key technical idea of ​​this invention is to reduce the dimensionality of the high-dimensional geometric parameter space to a sensitive parameter subspace through principal component analysis. This method identifies the parameter directions that have the most significant impact on aerodynamic performance based on the eigenvalue decomposition of the covariance matrix. While retaining the main optimization information, it compresses the search dimension from 20 to 30 dimensions to 3 to 5 dimensions. Compared with traditional full-parameter optimization methods, this significantly reduces the number of computational fluid dynamics simulations and avoids optimization search noise introduced by insensitive parameters, enabling the optimization algorithm to converge to the global optimum more efficiently. The second key technical idea is to construct an unsteady correction factor model to compensate for the prediction bias of steady-state calculations on the separated flow. This model maps the separation region feature parameters extracted by steady-state calculations to correction coefficients, obtaining a performance evaluation close to the unsteady real flow field while maintaining the computational efficiency of steady Reynolds-averaged Navier-Stokes equations. Compared with the method of directly using unsteady calculations, it saves tens of times the computation time, and compared with the traditional method that relies solely on steady-state calculations, it significantly improves the accuracy of the pressure recovery coefficient prediction in the diffuser section. The synergy between these two technical approaches lies in the fact that dimensionality reduction optimization ensures that the parameter space can be fully explored under limited computing resources, while unsteady correction guarantees the accuracy of each evaluation. The combination of the two enables the entire optimization process to achieve the optimal balance between computational efficiency and result reliability, ultimately realizing efficient and accurate optimization of the aerodynamic profile of the full-size subsonic wind tunnel loop.

[0063] It should be noted that this invention also solves the following technical problem: the high-dimensional geometric parameter space in the optimization of the aerodynamic profile of a full-size subsonic wind tunnel loop leads to a huge demand for computational fluid dynamics samples and low optimization efficiency. Traditional methods directly perform optimization searches within a complete parameter space consisting of twenty to thirty geometric parameters, including the length of the diffuser section, the outlet area of ​​the diffuser section, the corner radius, and the coordinates of the control points of the contraction section profile. The surrogate model requires a large number of sample points to establish a reliable response surface, and each sample point requires a time-consuming three-dimensional computational fluid dynamics solution, resulting in an unacceptable optimization cycle. This invention identifies the eigenvalues ​​of the covariance matrix of geometric parameters through principal component analysis, selects principal component directions with a cumulative contribution rate of eigenvalues ​​exceeding 0.85 to construct a sensitive parameter subspace, and compresses the original twenty to thirty-dimensional parameter space to three to five principal component dimensions that significantly affect the optimization objective. While maintaining optimization accuracy, the search dimension is reduced by more than 80%, reducing the number of initial sample points required by the Kriging surrogate model from hundreds to dozens, significantly reducing the number of computational fluid dynamics calls and improving optimization iteration efficiency, ensuring that the complex multi-objective optimization of the full-size wind tunnel loop is completed within an engineering-acceptable timeframe.

[0064] Specifically, the principle of this invention is as follows: The fundamental reason why this invention can solve the prediction bias problem of the steady-state calculation framework lies in establishing a bridging relationship between the flow separation characteristics obtainable by steady-state calculation and the unsteady real performance. Although the steady-state Reynolds-averaged Navier-Stokes equations cannot accurately predict the dynamic evolution of the separation zone, they can still provide predictions with relative trend significance for the geometric scale characteristics of the separation zone, such as the length of the separation zone and the height of the recirculation zone. These geometric characteristics have a definite physical correlation with the separation Reynolds number and the adverse pressure gradient parameter. The unsteady correction factor model establishes a functional mapping from the geometric characteristics to the correction coefficients through neural networks or polynomial regression fitting. The correction coefficients essentially quantify the systematic deviation of steady-state calculations relative to the unsteady real flow field in pressure recovery prediction. Applying the correction coefficients to the pressure recovery coefficient of the diffusion section in steady-state calculations can restore the performance indicators under real unsteady conditions. This allows the optimization algorithm to approximate the Pareto front of the real flow field rather than the pseudo-front generated by steady-state calculations when searching based on the corrected objective function, ensuring the engineering effectiveness of the optimization results.

[0065] The following provides a specific embodiment 1 of the present invention, and the specific implementation of each step in this embodiment 1 is described in detail below.

[0066] The specific implementation of step S01 involves acquiring the initial geometric parameters of the full-size subsonic wind tunnel loop, including the length of the diffuser section. The unit is meters (m), representing the outlet area of ​​the diffusion section. The unit is Corner radius The unit is meters (m), and the coordinates of the control points for the contraction section profile are as follows. The units are all in meters, where , The number of control points is dimensionless and typically ranges from 6 to 8; boundary condition parameters include the total inlet pressure. The unit is Pa, and the inlet total temperature is... Unit is K, outlet static pressure The unit is Pa. The above geometric parameters are obtained by measuring the actual dimensions of the wind tunnel loop, and the boundary condition parameters are determined according to the wind tunnel design conditions.

[0067] The specific implementation of step S02 involves performing principal component analysis on the initial geometric parameters, first constructing the initial geometric parameter vector. ,in This is a reference length, in meters, with a value of 1 meter. For reference area, the unit is... The value is 1 By normalizing the geometric parameter vector to make each component dimensionless, the covariance matrix is ​​calculated. The formula is expressed as follows:

[0068] ;

[0069] In the formula, The covariance matrix is ​​the first Line number Column elements, dimensionless; The sample size is dimensionless. For the first The first sample A normalized geometric parameter, dimensionless; For the first The mean of a normalized geometric parameter, dimensionless; For the first The first sample A normalized geometric parameter, dimensionless; For the first Find the mean of the normalized geometric parameters, which are dimensionless. Solve for the eigenvalues ​​of the covariance matrix. With feature vectors eigenvalues Dimensionless, eigenvector Each component is dimensionless, and the formula for the contribution rate of the eigenvalues ​​is as follows:

[0070] ;

[0071] In the formula, For the first The contribution rate of each principal component is dimensionless. For the first 1 eigenvalue, dimensionless; The original parameter dimension is dimensionless; The sum of all eigenvalues ​​is dimensionless. Identify those that satisfy... The principal component directions project the initial geometric parameters onto the sensitive parameter subspace, and the formula for the dimensionality-reduced parameter vector is as follows:

[0072] ;

[0073] In the formula, The parameter vector is a dimensionless vector. For the reason before The projection matrix consists of eigenvectors, where each element is dimensionless. The number of principal components to be retained is dimensionless and typically ranges from 3 to 5. This is the initial normalized geometric parameter vector, and each component is dimensionless.

[0074] The specific implementation of step S03 involves generating an initial sample point set using Latin hypercube sampling within the dimension-reduced parameter vector space, and defining the value range of each dimension-reducing parameter. Divided into equal parts There are several intervals, ensuring that one sample point is drawn from each interval. For the first The minimum value of each dimension reduction parameter, dimensionless. For the first The maximum value of each dimension reduction parameter, dimensionless. The total number of sample points is dimensionless and typically ranges from 30 to 50. A three-dimensional computational domain mesh is constructed for each sample point in the initial sample point set, and the steady Reynolds-averaged Navier-Stokes equations are solved to extract the uniformity of the exit velocity in the contraction section. pressure recovery coefficient of the diffusion section Total pressure loss coefficient in corner area As a response value, the formula for the uniformity of the exit velocity of the contraction section is expressed as follows:

[0075] ;

[0076] In the formula, The uniformity of the exit velocity of the contraction section is dimensionless. The standard deviation of the velocity at the exit section of the contraction section is expressed in m / s. The average velocity at the exit section of the contraction section is expressed in m / s. The reference speed standard deviation is expressed in m / s, with an empirical value of 5 m / s. The reference average velocity is expressed in m / s, with an empirical value of 50 m / s. The formula for the pressure recovery coefficient in the diffuser section is as follows:

[0077] ;

[0078] In the formula, The pressure recovery coefficient of the diffusion section is dimensionless. This is the static pressure at the outlet of the diffuser section, in Pa. This refers to the static pressure at the inlet of the diffuser section, expressed in Pa. For reference air density, the unit is kg / The value is 1.225 kg / ; The reference speed is in m / s, and the value is taken as 50 m / s. The formula for the total pressure loss coefficient in the corner area is as follows:

[0079] ;

[0080] In the formula, The total pressure loss coefficient in the corner area is dimensionless. The total pressure at the inlet of the corner area is expressed in Pa. The total pressure at the outlet of the corner area is expressed in Pa. The reference speed for the corner area is in m / s, with an empirical value of 50 m / s.

[0081] The specific implementation of step S04 involves constructing a Kriging surrogate model using the dimensionality-reduced parameter vector and response values ​​of the initial sample point set. The formula for the predicted value of the Kriging model is expressed as follows:

[0082] ;

[0083] In the formula, Design point The Kriging prediction value at the point can be any response value among the uniformity of the outlet velocity of the contraction section, the pressure recovery coefficient of the diffusion section, or the total pressure loss coefficient of the corner region, and is dimensionless. The mean of the response values ​​is dimensionless. Let be the correlation vector between the point to be predicted and the known sample points. Each component is dimensionless, and its i-th... The formula for calculating each element is: ,in For the first The relevant parameters are dimensionless and determined through maximum likelihood estimation. This is a smoothing parameter, dimensionless, and typically takes a value of 2. The first point to be predicted Each dimension-reduced parameter component is dimensionless. For the first The first sample point Each dimension-reduced parameter component is dimensionless. Given the correlation matrix between sample points, where each element is dimensionless, its i.e., the nth element... Line number The formula for calculating column elements is: ,in and These are the row and column indices of the correlation matrix, respectively; they are dimensionless. For the first The first sample point Each dimension-reduced parameter component is dimensionless. For the first The first sample point Each dimension-reduced parameter component is dimensionless. The response value vector of known sample points, with each component being dimensionless; Given a vector of all 1s with dimensionless components, calculate the root mean square error between the Kriging surrogate model predictions and the computational fluid dynamics calculations. The formula for the root mean square error is as follows:

[0084] ;

[0085] In the formula, The root mean square error is dimensionless. The number of test samples is dimensionless and typically ranges from 10 to 20. For the first Kriging predictions for each test point, dimensionless; For the first The computational fluid dynamics values ​​for each test point are dimensionless. This is the average value of the test sample response, which is dimensionless. If the root mean square error is less than 5%, proceed to step S05; otherwise, add supplementary sample points, update the Kriging surrogate model, and repeat step S04.

[0086] The specific implementation of step S05 is to calculate the desired improved function value based on the Kriging surrogate model. The formula for the desired improved function is expressed as follows:

[0087] ;

[0088] In the formula, Design point The expected improved function value at point is dimensionless; The current optimal objective function value is dimensionless. These are Kriging predictions, dimensionless. The standard deviation of the Kriging prediction is dimensionless and is calculated using the following formula: ,in The standard deviation of the response values ​​is dimensionless. The cumulative distribution function is a standard normal distribution; the input parameters are dimensionless, and the output value is dimensionless. Given a standard normal probability density function, the input parameters and output values ​​are dimensionless. The top 20 candidate points with the largest expected improvement function values ​​are selected as the initial positions for a multi-starting point parallel search. For each initial position, the steady Reynolds-averaged Navier-Stokes equation is solved, and the separation zone length is extracted. Height of the reflux zone ,in The unit is meters. The unit is meters (m).

[0089] The specific implementation of step S06 involves inputting the separation zone length and reflux zone height into the unsteady correction factor model, and calculating the separation Reynolds number using the unsteady correction factor model. With inverse pressure gradient parameters The formula for separating the Reynolds number is expressed as follows:

[0090] ;

[0091] In the formula, To separate the Reynolds number, dimensionless; The boundary layer momentum thickness at the separation point, in meters, is obtained through boundary layer integration. The calculation formula is as follows: ,in The boundary layer thickness, in meters, is determined by the location where the velocity distribution reaches 99% of the incoming flow velocity. Distance from the wall The velocity at that point is expressed in m / s. This is the distance from the wall, in meters (m). The velocity of the incoming flow is expressed in m / s. Kinematic viscosity, unit: / s, for air under standard conditions, takes the value of / s. / s. The formula for the inverse pressure gradient parameter is as follows:

[0092] ;

[0093] In the formula, The inverse pressure gradient parameter is dimensionless. This represents the static pressure difference along the friction path, in Pa. The characteristic length is measured in meters (m) and is typically taken as the length of the diffusion segment. The reference length is in meters (m), and its value is 1 meter. The formula for the correction factor is as follows:

[0094] ;

[0095] In the formula, This is a correction factor, dimensionless; The weights for separating the Reynolds number are dimensionless and have an empirical value of 0.4. The weights for correcting the inverse pressure gradient parameters are dimensionless and have an empirical value of 0.6. For reference separation, the Reynolds number is dimensionless and typically takes the value of . ; The reference adverse pressure gradient parameter is dimensionless and typically takes a value of 0.02. The formula for the corrected pressure recovery coefficient in the diffuser section is as follows:

[0096] ;

[0097] In the formula, The corrected pressure recovery coefficient for the diffusion section is dimensionless. The pressure recovery coefficient of the diffusion section is obtained by solving the steady Reynolds-averaged Navier-Stokes equation in step S05, and is dimensionless.

[0098] The specific implementation of step S07 is to obtain the corrected diffusion section pressure recovery coefficient. The uniformity of the exit velocity of the contraction section extracted in step S05 The total pressure loss coefficient of the corner area extracted in step S05 The three performance indicators mentioned above are normalized by dividing them by their corresponding reference values. The formula for the normalized uniformity of the exit velocity of the contraction section is expressed as follows:

[0099] ;

[0100] In the formula, The normalized uniformity of the exit velocity of the contraction section is dimensionless. For reference purposes regarding the uniformity of the exit velocity in the contraction section, a dimensionless value is typically taken as 0.95. The formula for the normalized corrected pressure recovery coefficient in the diffusion section is as follows:

[0101] ;

[0102] In the formula, The normalized pressure recovery coefficient of the diffusion section is dimensionless. The pressure recovery coefficient for the reference diffusion section is dimensionless and typically taken as 0.75. The formula for the normalized total pressure loss coefficient in the corner region is as follows:

[0103] ;

[0104] In the formula, The normalized total pressure loss coefficient for the corner region is dimensionless. The total pressure loss coefficient in the reference corner area is dimensionless and determined based on wind tunnel design experience, with an empirical value of 0.15. The formula for the multi-objective optimization function is as follows:

[0105] ;

[0106] In the formula, The function is a multi-objective optimization function, dimensionless; the weight coefficients 0.35, 0.45, and 0.20 are calculated using the analytic hierarchy process (AHP) and are also dimensionless. A non-dominated sorting genetic algorithm is used to search for the Pareto front solution set within the dimension-reduced parameter vector space.

[0107] The specific implementation of step S08 is to calculate the comprehensive performance index of each solution in the Pareto front solution set. The formula for the comprehensive performance index is as follows:

[0108] ;

[0109] In the formula, This is a comprehensive performance index, dimensionless. The parameter combination with the highest comprehensive performance index is selected. The parameter combination is inversely transformed from the sensitive parameter subspace to the complete geometric parameter space. The formula for the inverse transformation is as follows:

[0110] ;

[0111] In the formula, The optimized parameter vector is a complete geometric parameter space, and each component is dimensionless. This is a projection matrix, and all its elements are dimensionless. The optimized, dimensionless parameter vector has dimensionless components. A hierarchical parameterization method is used to reconstruct the aerodynamic profile of the wind tunnel loop.

[0112] The specific implementation of step S09 involves performing steady Reynolds-averaged Navier-Stokes equations to verify the reconstructed wind tunnel loop aerodynamic profile, and obtaining the uniformity of the exit velocity in the contraction section from the verification calculation. With the corrected pressure recovery coefficient of the diffuser section ,in Dimensionless Dimensionless, when the calculated uniformity of the exit velocity in the contraction section is greater than 0.95 and the calculated pressure recovery coefficient in the corrected diffusion section is greater than 0.75, the reconstructed wind tunnel loop aerodynamic profile is output as the optimization result; otherwise, the separation Reynolds number correction weight in the unsteady correction factor model is adjusted. Weights adjusted with inverse pressure gradient parameters Then return to step S06, and adjust the strategy so that when the verification result does not meet the threshold, the separation Reynolds number correction weight is increased by 0.05 to 0.1, and the inverse pressure gradient parameter correction weight is reduced accordingly to keep the weight sum at 1. The adjustment formula is expressed as follows:

[0113] ;

[0114] ;

[0115] In the formula, The weights are adjusted for the separated Reynolds number; they are dimensionless. The weights are adjusted for the inverse pressure gradient parameters and are dimensionless. This is the weight adjustment amount, dimensionless, with a value range of 0.05 to 0.1.

[0116] The specific implementation of principal component analysis is to construct the covariance matrix of the initial geometric parameters and solve for the eigenvalues ​​and eigenvectors. The formula for the covariance matrix is ​​the same as in step S02, and will not be elaborated here. When the ratio of the eigenvalue of a principal component to the sum of all eigenvalues ​​exceeds 0.85, the principal component is considered to contain the main information of the original parameters, thus compressing the high-dimensional space containing dozens of parameters into a low-dimensional sensitive parameter subspace composed of 3 to 5 principal components.

[0117] The specific implementation method for the sensitive parameter subspace is the same as in step S02, and will not be described in detail here. The specific implementation method for Latin hypercube sampling is the same as in step S03, and will not be described in detail here. The specific implementation method for the Kriging surrogate model is the same as in step S04, and will not be described in detail here. The specific implementation method for improving the function value is the same as in step S05, and will not be described in detail here. The specific implementation method for multi-starting point parallel search is the same as in step S05, and will not be described in detail here.

[0118] The specific implementation of the unsteady correction factor model involves establishing empirical correction relationships through calculations using a small number of unsteady Reynolds-averaged Navier-Stokes equations or separated eddy simulations. This process addresses the separation zone length, which cannot be accurately predicted in steady calculations. Height of the reflux zone As input variables, the functional relationship between the correction coefficient and the separation zone length and the recirculation zone height is fitted by neural networks or polynomial regression. The correction coefficient is then multiplied by the pressure recovery coefficient of the diffusion section calculated in steady condition to compensate for the systematic bias caused by the insufficient prediction of the separation flow by the steady Reynolds-averaged Navier-Stokes equations.

[0119] The specific implementation method for separating the Reynolds number is the same as in step S06, and will not be described in detail here, wherein the boundary layer momentum thickness The calculation formula is The specific implementation of the inverse pressure gradient parameters is the same as in step S06, and will not be described in detail here. The specific implementation of the non-dominated sorting genetic algorithm is the same as in step S07, and will not be described in detail here. The specific implementation of the Pareto front solution set is the same as in step S07, and will not be described in detail here.

[0120] The specific implementation of the Analytic Hierarchy Process (AHP) involves constructing a judgment matrix to express the relative importance of different objectives. The formula is expressed as follows:

[0121] ;

[0122] In the formula, the judgment matrix Each element is dimensionless. The element in the first row and second column represents the importance ratio of the outlet velocity uniformity of the contraction section to the pressure recovery coefficient of the diffusion section, which is 1:1.3. The element in the first row and third column represents the importance ratio of the outlet velocity uniformity of the contraction section to the total pressure loss coefficient in the corner region, which is 1.75:1. The element in the second row and third column represents the importance ratio of the pressure recovery coefficient of the diffusion section to the total pressure loss coefficient in the corner region, which is 2.25:1. Calculate the largest eigenvalue of the judgment matrix. Corresponding feature vector And normalize, where Dimensionless Each component is dimensionless, and the formula for the normalized weight vector is as follows:

[0123] ;

[0124] In the formula, For the first The normalized weighting coefficients of each objective, dimensionless; For feature vectors The One component, dimensionless; The sum of all components of the eigenvector is dimensionless. The weight coefficients obtained from the above calculations are 0.35, 0.45, and 0.20.

[0125] The specific implementation of the hierarchical parameterization method is to reconstruct the wind tunnel loop profile using a multi-scale geometric expression strategy. The large-scale profile defines the overall direction of the diffusion and contraction sections through 6 to 8 B-spline curve control points. Small-scale local features are superimposed on the free-deformable mesh to apply local curvature adjustment in the corner transition area and the starting section of the diffusion section. The displacement amplitude of the control points is limited to ±15% of the initial position to prevent unreasonable concavity and convexity abrupt changes.

[0126] The specific implementation methods for separating the Reynolds number correction weight and the adverse pressure gradient parameter correction weight are the same as in steps S06 and S09, and will not be described in detail here. The threshold value of 0.95 for the uniformity of the exit velocity in the contraction section is derived from the engineering requirements for the flow field quality in the full-scale subsonic wind tunnel test section. The threshold value of 0.75 for the pressure recovery coefficient in the diffusion section is derived from the wind tunnel loop energy loss control standard.

[0127] To better understand and implement this invention, the following is a specific application scenario of this invention, Example 2:

[0128] The test section of a wind tunnel has a cross-sectional dimension of 8×6m and a design wind speed range of 20 to 100m / s. The original loop had problems such as a low pressure recovery coefficient in the diffuser section and uneven velocity distribution at the exit of the contraction section, resulting in the drive power exceeding the design value by 18% and the flow field quality in the test section failing to meet the requirements for precision measurement. The technical team decided to redesign the aerodynamic profile of the loop using a multi-objective optimization method based on steady CFD.

[0129] The technical team first collected the initial geometric parameters of the existing wind tunnel loop, including the lengths of the four diffuser sections of 12.5m, 15.8m, 14.2m, and 16.3m, corresponding to exit areas of 78m and 78m respectively. 95 112 128 The four corner radii are 5.2m, 4.8m, 5.5m, and 4.9m, respectively. The contraction section profile is defined by 18 control points, with coordinates ranging from 0 to 22m in the x-direction and from -4 to 4m in the y-direction. Boundary conditions are set as follows: inlet total pressure 101325Pa, inlet total temperature 288K, and outlet static pressure adjusted to the corresponding static pressure value according to wind speed requirements. The technical team performed principal component analysis on these 28 initial geometric parameters, constructed the covariance matrix, and solved for the eigenvalues. They identified that the ratios of the eigenvalues ​​of the first four principal components to the sum were 0.42, 0.28, 0.19, and 0.11, respectively, with a cumulative contribution rate of 0.89, exceeding the threshold of 0.85. The technical team projected the original 28-dimensional parameters onto the sensitive parameter subspace formed by the directions of these four principal components, forming a dimension-reduced parameter vector. .

[0130] Within the reduced-dimensional parameter vector space, the technical team used Latin hypercube sampling to generate 65 initial sample points, with each parameter dimension equally divided into 65 intervals to ensure uniform sample distribution. For each sample point, the team constructed a 3D computational domain containing 3.4 million hexahedral mesh elements, with the mesh refined near the walls to a first-layer mesh height of 0.08mm to meet certain requirements. The requirement is less than 1. Solve the steady Reynolds-averaged Navier-Stokes equations using Realizable. The turbulence model sets the convergence criterion as the residuals of each equation being less than 1. The technical team extracted the uniformity of the exit velocity in the contraction section, the pressure recovery coefficient in the diffusion section, and the total pressure loss coefficient in the corner region from the calculation results as response values, as shown in Table 1.

[0131] Table 1 Aerodynamic performance response values ​​of some sample points

[0132]

[0133] The technical team constructed a Kriging surrogate model using the dimensionality-reduced parameter vector and response values ​​of 65 sample points, selected the Gaussian correlation function as the spatial covariance structure, and determined the correlation length parameter through maximum likelihood estimation. The technical team calculated the root mean square error between the surrogate model's predicted values ​​and the CFD calculation values. The error in the uniformity of the velocity at the exit of the contraction section was 3.2%, the error in the pressure recovery coefficient of the diffusion section was 4.1%, and the error in the total pressure loss coefficient of the corner region was 3.8%, all of which were less than the 5% threshold, indicating that the surrogate model has sufficient prediction accuracy.

[0134] The technical team calculated the expected improvement function value based on the Kriging surrogate model, generated 50,000 candidate points through Monte Carlo sampling, evaluated their expected improvement values, and selected the top 20 candidate points with the largest expected improvement function values ​​as the initial positions for multi-starting point parallel search. For these 20 initial positions, the team performed parallel solutions to the steady Reynolds-averaged Navier-Stokes equations, identified the separation zone within the diffusion section from the flow field results, and extracted the length of the separation zone and the height of the recirculation zone, such as... Figure 2 As shown. The typical separation zone length ranges from 1.8 to 4.6 m, and the recirculation zone height ranges from 0.15 to 0.52 m.

[0135] The technical team input the separation zone length and recirculation zone height into a pre-established unsteady correction factor model. This model, trained via a neural network based on 35 sets of unsteady Reynolds-averaged Navier-Stokes equations, contains 8 hidden neurons. The model first calculates the separation Reynolds number, taking the boundary layer momentum thickness at the separation point as 0.024 m, the incoming flow velocity as 85 m / s, and the kinematic viscosity as... The separated Reynolds number was obtained as 136000. Subsequently, the adverse pressure gradient parameters were calculated, taking a static pressure difference of 2850 Pa, a dynamic pressure of 3612 Pa, and a characteristic length of 14.2 m, resulting in an adverse pressure gradient parameter of 0.056. Based on the separated Reynolds number and the adverse pressure gradient parameters, the model output correction coefficients ranged from 1.12 to 1.18. The technical team multiplied these correction coefficients by the steady-state calculated pressure recovery coefficient of the diffusion section to obtain the corrected pressure recovery coefficient of the diffusion section, as shown below. Figure 3 As shown, the corrected values ​​are on average 14% higher than the steady-state calculation results.

[0136] After obtaining the corrected pressure recovery system in the diffuser section, the uniformity of the outlet velocity in the contraction section, and the total pressure loss coefficient in the corner region, the technical team normalized the three objectives by dividing them by reference values ​​using weighting coefficients determined by the analytic hierarchy process (AHP). The reference values ​​for the uniformity of the outlet velocity in the contraction section were set to 0.95, the pressure recovery coefficient in the diffuser section to 0.78, and the total pressure loss coefficient in the corner region to 0.045. A multi-objective optimization function was then constructed. A non-dominated sorting genetic algorithm was used to search for the Pareto front solution set in a dimension-reduced parameter vector space. The technical team set the population size to 120, the number of generations to 80, the crossover probability to 0.85, and the mutation probability to 0.12. After iterative calculation, a Pareto front containing 58 non-dominated solutions was obtained.

[0137] The technical team calculated the comprehensive performance index for each solution in the Pareto front by multiplying the normalized exit velocity uniformity of the contraction section by 0.35, the corrected pressure recovery coefficient of the diffusion section by 0.45, and the reciprocal of the total pressure loss coefficient in the corner region by 0.20, and then summing the results. The parameter combination with the highest comprehensive performance index corresponds to... , , , The technical team inversely transformed the parameter combination from the sensitive parameter subspace to the complete geometric parameter space, obtaining optimized diffuser lengths of 13.8m, 16.5m, 15.1m, and 17.2m, and corner radii of 5.8m, 5.2m, 6.1m, and 5.5m, respectively. The coordinates of the control points for the contraction section profile were adjusted accordingly. Using a hierarchical parameterization method, the team defined the large-scale profile of the diffuser section using seven B-spline control points. Local curvature adjustments were applied by overlaying a free-deformable mesh in the corner transition zone, with the curvature radius continuously varying from 3.2 to 6.8m. The control point displacement was limited to ±15% of the initial position, reconstructing the complete aerodynamic profile of the wind tunnel loop.

[0138] The technical team performed steady Reynolds-averaged Navier-Stokes equations to verify the reconstructed wind tunnel loop aerodynamic profile, increasing the mesh size to 5.2 million to improve analytical accuracy. The verification calculations showed an exit velocity uniformity of 0.962 in the contraction section, exceeding the set threshold of 95%, and a pressure recovery coefficient of 0.782 in the corrected diffuser section, exceeding the set threshold of 0.75. Both indicators meet engineering requirements. The technical team output the reconstructed wind tunnel loop aerodynamic profile as the final optimization result.

[0139] The advancements of this invention over traditional methods are reflected in three aspects. First, by using principal component analysis, the high-dimensional geometric parameter space is compressed into a sensitive parameter subspace, avoiding the waste of computational resources caused by a large number of insensitive parameters in traditional full-parameter space optimization. This allows the optimization search to focus on parameter directions that significantly affect the objective function, significantly reducing the number of CFD calculations required. Second, the introduction of an unsteady correction factor model compensates for the systematic bias of the steady Reynolds-averaged Navier-Stokes equations in predicting separated flows. While maintaining the efficiency of steady calculations, it obtains performance evaluations close to the unsteady real flow field, solving the problem of unreliable optimization results due to the overestimation of the pressure recovery coefficient in the diffusion section under the traditional steady calculation framework. Third, the sequential addition strategy using a Kriging surrogate model combined with the expected improvement criterion balances the needs of global exploration and local development. Compared with the traditional genetic algorithm that directly calls CFD calculations, this significantly reduces the risk of getting trapped in local optima, ensuring that the search process can approach the global Pareto front.

[0140] It should be noted that the variables involved in this invention are explained in detail in Tables 2 and 3.

[0141] Table 2. Variable Explanation Table (Part 1)

[0142]

[0143] Table 3. Variable Explanation Table (Part Two)

[0144]

[0145] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in the present invention should be included within the scope of protection of the present invention.

Claims

1. A multi-objective optimization method for the aerodynamic profile of a full-size subsonic wind tunnel loop based on steady CFD, characterized in that, include: Collect initial geometric parameters and set boundary condition parameters; Principal component analysis is performed on the initial geometric parameters to form a dimension-reduced parameter vector; Latin hypercube sampling is used to generate an initial sample point set in the dimension-reduced parameter vector space, and the steady Reynolds-averaged Navier-Stokes equation is executed to extract the response value; a Kriging surrogate model is constructed; the expected improvement function value is calculated based on the Kriging surrogate model, and a multi-starting-point parallel search is performed to extract the separation zone length and recirculation zone height; The separation zone length and recirculation zone height are input into the unsteady correction factor model to calculate the correction coefficient and obtain the pressure recovery coefficient of the diffusion section after correction; a multi-objective optimization function is constructed and a non-dominated sorting genetic algorithm is used to search for the Pareto front solution set; the comprehensive performance index is calculated, the optimal parameter combination is selected, and the aerodynamic profile of the wind tunnel loop is reconstructed; the verification calculation is performed and the optimization results are output.

2. The method according to claim 1, characterized in that, The initial geometric parameters are specifically the length of the diffuser section, the outlet area of ​​the diffuser section, the corner radius, and the coordinates of the control points of the contraction section profile. The boundary condition parameters are specifically the total inlet pressure, the total inlet temperature, and the static outlet pressure.

3. The method according to claim 2, characterized in that, The steps of principal component analysis specifically involve constructing the covariance matrix of the initial geometric parameters and solving for eigenvalues ​​and eigenvectors to achieve dimensionality reduction, identifying principal component directions where the ratio of eigenvalue to the sum of eigenvalues ​​is greater than 0.85, and projecting the initial geometric parameters onto the sensitive parameter subspace.

4. The method according to claim 3, characterized in that, The sensitive parameter subspace refers to the low-dimensional parameter space that has a significant impact on the optimization objective, which is screened out through principal component analysis, reducing the search dimension of the optimization problem from the original 20 to 30 dimensions to 3 to 5 dimensions.

5. The method according to claim 4, characterized in that, The Latin hypercube sampling refers to a stratified sampling technique that divides the range of values ​​for each parameter into several intervals and extracts exactly one sample from each interval.

6. The method according to claim 5, characterized in that, The step of extracting response values ​​by solving the steady Reynolds-averaged Navier-Stokes equations involves constructing a three-dimensional computational domain grid for each sample point in the initial sample point set and solving the steady Reynolds-averaged Navier-Stokes equations to extract the uniformity of the exit velocity in the contraction section, the pressure recovery coefficient in the diffusion section, and the total pressure loss coefficient in the corner region as response values.

7. The method according to claim 6, characterized in that, The Kriging agent model refers to an alternative model based on spatial interpolation theory, which assumes that the response value is composed of a deterministic trend function and a stochastic process, and determines the hyperparameters of the correlation function through maximum likelihood estimation.

8. The method according to claim 7, characterized in that, The steps for constructing the Kriging surrogate model are as follows: the Kriging surrogate model is constructed using the dimensionality-reduced parameter vector and response values ​​of the initial sample point set; the root mean square error between the predicted value of the Kriging surrogate model and the calculated value of the computational fluid dynamics is calculated; if the root mean square error is less than 5%, proceed to the next step; otherwise, supplementary sample points are added and the Kriging surrogate model is updated.

9. The method according to claim 8, characterized in that, Specifically, the expected improvement function value is determined by comprehensively considering the gap between the Kriging surrogate model prediction value and the current optimal objective function value, as well as the prediction uncertainty, and selecting the top 20 candidate points with the largest expected improvement function values ​​as the initial positions for multi-starting point parallel search.

10. The method according to claim 9, characterized in that, The step of performing a multi-starting-point parallel search to extract the separation zone length and the recirculation zone height specifically involves solving the steady Reynolds-averaged Navier-Stokes equations at each initial position and extracting the separation zone length and the recirculation zone height. The separation zone length and the recirculation zone height are used as flow separation characteristic quantities under the steady calculation framework.

Images