A method for aerodynamic optimization design of wind tunnel internal surfaces

Abstract

This invention provides a method for aerodynamic optimization design of wind tunnel profiles, belonging to the field of wind tunnel profile technology. This invention identifies high-gradient feature points by extracting the curvature distribution of profile features, generates a hybrid mesh using a graph-driven adaptive topology reconstruction algorithm, performs flow field calculations based on the hybrid RANS-LES turbulence simulation method, solves the control equations using a high-precision shock capture format coupled with a WENO limiter, and performs adaptive mesh refinement. It adjusts the coordinates of profile control points using a topological potential-based asymptotic optimization algorithm, and determines whether to output the optimization results or return to iterate again by judging whether the aerodynamic performance evaluation parameters meet the convergence criteria. This solves the technical problem of low optimization efficiency and difficulty in converging to the optimal solution caused by the disconnect between mesh generation, flow field calculation, and profile adjustment in the aerodynamic optimization process of wind tunnel profiles.

Description

Technical Field

[0001] This invention belongs to the field of wind tunnel internal surface technology, and more specifically, relates to a method for aerodynamic optimization design of wind tunnel internal surfaces. Background Technology

[0002] Wind tunnel aerodynamic design of airfoils is a key technology in the aerospace field. Traditional methods determine the airfoil shape through empirical formulas and simplified flow theories, and then use computational fluid dynamics numerical simulations to verify aerodynamic performance. In current wind tunnel design practice, designers typically first establish an initial airfoil geometry model, then use commercial software to generate a computational mesh, then call a flow field solver for numerical simulation, and finally manually adjust the airfoil parameters based on the calculation results. This process requires multiple iterations to obtain a satisfactory design. Because the mesh generation uses a fixed topology, it cannot adapt to changes in flow field characteristics; the flow field calculation using a single turbulence model struggles to balance computational accuracy and efficiency; and airfoil optimization relies on the designer's experience and lacks a systematic adjustment strategy, resulting in each stage of the optimization process being independent and unable to form a closed-loop feedback mechanism. In other words, existing technologies suffer from the technical problem of low optimization efficiency and difficulty in converging to the optimal solution due to the disconnect between mesh generation, flow field calculation, and airfoil adjustment in the wind tunnel aerodynamic optimization process. Summary of the Invention

[0003] In view of this, the present invention provides a method for optimizing the aerodynamic design of a wind tunnel profile, which can solve the technical problem in the prior art where the three stages of mesh generation, flow field calculation and profile adjustment are isolated from each other, resulting in low optimization efficiency and difficulty in converging to the optimal solution.

[0004] This invention is implemented as follows: It provides a method for aerodynamic optimization design of wind tunnel surfaces. An initial geometric model of the wind tunnel surface is established, and the surface curvature distribution is extracted. High-gradient feature points are marked in regions with drastic curvature changes. A graph-driven adaptive topology reconstruction algorithm is used to generate a hybrid mesh, generating a prism layer mesh near the wall and a polyhedral mesh in the core region. Numerical flow field calculations are performed based on a hybrid RANS-LES turbulence simulation method. Large eddy simulation is used in the flow separation prediction region, and the Reynolds average method is used in the attached flow region. A high-precision shock wave capture scheme coupled with a WENO limiter is used to solve the governing equations. Adaptive mesh refinement is performed in the shock wave and boundary layer interference region. Pressure and velocity distribution data of the surface are extracted as aerodynamic performance evaluation parameters. A surface asymptotic optimization algorithm based on topological potential is used to adjust the coordinates of the surface control points. The aerodynamic performance evaluation parameters are judged to meet the convergence criteria. If they do, the optimized wind tunnel surface geometric model is output; otherwise, the mesh generation and flow field calculation are repeated until convergence.

[0005] The identification criterion for high-gradient feature points is that the difference in surface curvature between adjacent control points is greater than 0.15. .

[0006] Specifically, the graph theory-driven adaptive topology reconstruction algorithm abstracts the computational grid into a graph structure, where grid nodes correspond to graph vertices and cell connection relationships correspond to graph edges, and identifies regions with poor grid quality through graph connectivity analysis.

[0007] In the graph theory-driven adaptive topology reconstruction algorithm, the graph cut algorithm is used to separate the sub-regions that need to be reconstructed, the inter-layer connectivity of the prism layer grid is optimized based on the shortest path algorithm, the minimum spanning tree algorithm is used to determine the priority order of grid densification, and the spectral clustering method is introduced to aggregate units with similar flow characteristics.

[0008] The graph theory-driven adaptive topology reconstruction algorithm employs a three-tiered precision approach for progressive optimization. The first tier of precision corresponds to the ratio of the grid cell size to the wind tunnel feature length in the coarse grid stage. ∈[0.05, 0.10], the second-order precision corresponds to the ratio of the grid cell size to the wind tunnel feature length in the medium-grid stage. ∈[0.01, 0.05), the third-order precision corresponds to the ratio of the grid cell size to the wind tunnel characteristic length in the fine-grid stage. ∈[0.001, 0.01).

[0009] The height of the first layer of the prism mesh is calculated based on the wall Reynolds number. The value is controlled within 1, and the number of prism mesh layers ranges from 15 to 30 layers.

[0010] In the hybrid RANS-LES turbulence simulation method, the Reynolds-averaged Navier-Stokes equations are used to solve the time-averaged flow field in the attached flow region. The Reynolds stress term in the Reynolds-averaged Navier-Stokes equations is obtained through SST. The turbulence model is closed, and large eddy simulation is used to directly solve the large-scale turbulence structure in the flow separation prediction region, while small-scale turbulence is modeled through a subgrid model.

[0011] The method for identifying the flow separation prediction region is to calculate the wall shear stress and pressure gradient parameters. When the wall shear stress is close to zero and the pressure gradient parameter is greater than the critical pressure gradient parameter threshold, it is determined to be a flow separation prediction region. The critical pressure gradient parameter threshold is 0.05.

[0012] The hybrid RANS-LES turbulence simulation method employs a two-step accuracy approach to optimize computational efficiency. In the first-step accuracy stage, the Reynolds average method is used globally to calculate the ratio of the initial solution time step to the flow characteristic time. ∈[10, 50], the second-order accuracy stage switches to the ratio of the large eddy simulation method time step to the flow characteristic time in the identified flow separation prediction region. ∈[0.1, 1.0].

[0013] Specifically, the high-precision shock wave capture format coupled with the WENO limiter uses the AUSM+ format to calculate the convection flux of the cell interface and a fifth-order WENO limiter to reconstruct the flow variables of the cell interface.

[0014] The criterion for adaptive mesh encryption is based on a density gradient sensor. When the density gradient sensor value is greater than 0.01, adaptive mesh encryption is triggered. The size of the mesh cells after adaptive mesh encryption is 0.5 times that of the original mesh cells.

[0015] The high-precision shock wave capture scheme coupled with the WENO limiter employs a three-step accuracy method to improve shock wave resolution. In the first-step accuracy stage, a second-order upwind scheme is used to roughly calculate the residual convergence to the initial residual. In the second-order accuracy stage, the AUSM+ scheme is switched to coupled third-order WENO limiter in the shock region, and the residual converges to the initial residual. In the third-order accuracy stage, a fifth-order WENO limiter is used in the shock wave and boundary layer interference region, and adaptive mesh refinement is performed to bring the residuals to converge to the initial residuals. times.

[0016] Specifically, the topological potential field-based asymptotic optimization algorithm for the shape is to construct a field theory model with pressure distribution as potential energy, regard each control point on the shape as a mass in the field, calculate the driving force of the potential energy field based on the local pressure gradient to drive the movement of the control point, introduce a topological sensing mechanism to identify flow separation regions and reattachment regions, and add virtual repulsive forces in the flow separation regions and reattachment regions.

[0017] The surface-based progressive optimization algorithm, based on topological potential field, employs a four-step accuracy approach to achieve progressive optimization from coarse to fine. The potential field intensity coefficient in the first-step accuracy stage... The ratio of the radius of action of the potential field of 1.0 to the total length of the profile. The ratio of the 0.20 control point movement step size to the surface feature size. The potential field intensity coefficient is 0.05, representing the second-order precision stage. The ratio of the radius of action of the potential field to the total length of the profile is 0.5. The ratio of the 0.10 control point movement step size to the surface feature size. The potential field strength coefficient is 0.02, representing the third-order precision stage. The ratio of the radius of action of the potential field to the total length of the profile is 0.2. The ratio of the 0.05 control point movement step size to the surface feature size. The potential field intensity coefficient is 0.005, representing the fourth-order precision stage. The ratio of the radius of action of the 0.05 potential field to the total length of the profile. The ratio of the 0.02 control point movement step size to the surface feature size. It is 0.001.

[0018] The aerodynamic performance evaluation parameters include pressure distribution uniformity index and velocity distribution stability index. The pressure distribution uniformity index is obtained by calculating the ratio of the standard deviation to the average value of the pressure coefficient of the profile, and the velocity distribution stability index is obtained by calculating the ratio of the maximum deviation of the axial velocity of the test section to the average velocity.

[0019] The convergence criterion is that the rate of change of aerodynamic performance evaluation parameters in three consecutive iterations is less than 0.5% and the maximum adjustment of the coordinates of the surface control points is less than 0.05% of the surface feature size.

[0020] This invention establishes a fully automated technical system encompassing feature recognition, mesh generation, flow field calculation, and profile optimization. It organically couples graph theory-driven adaptive topology reconstruction, hybrid turbulence simulation methods, and profile optimization algorithms based on topological potential fields. This achieves a collaborative optimization mechanism where the mesh dynamically adjusts according to flow field characteristics, computational accuracy progressively improves with each optimization stage, and profile adjustment automatically evolves along the pressure gradient direction. The invention employs a three- or four-step accuracy progression strategy, rapidly capturing large-scale flow characteristics in the coarse mesh stage, accurately locating shock waves and separation regions in the medium mesh stage, and finely distinguishing boundary layer interactions in the fine mesh stage. By combining dynamic control of mesh refinement, turbulence model switching, and profile adjustment step size, the optimization process naturally converges to the optimal solution that meets aerodynamic performance requirements. In summary, this invention solves the technical problem mentioned in the background art where the three stages of mesh generation, flow field calculation, and profile adjustment in wind tunnel profile aerodynamic optimization are isolated, leading to low optimization efficiency and difficulty in converging to the optimal solution. 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 surface pressure coefficient in the embodiment.

[0023] Figure 3 This is a distribution diagram of the density gradient sensors 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 of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings.

[0025] like Figure 1 The diagram shown is a flowchart of an aerodynamic optimization design method for wind tunnel surfaces provided by this invention. This method includes the following steps:

[0026] S01. Establish the initial geometric model of the wind tunnel surface and extract the surface feature curvature distribution. Mark high gradient feature points in areas with drastic curvature changes and output the coordinate set of high gradient feature points and surface feature curvature distribution data.

[0027] S02. A graph-driven adaptive topology reconstruction algorithm is used to generate a hybrid mesh. Prismatic layer meshes are generated near the wall, and polyhedral meshes are generated in the core region. The coordinates of the hybrid mesh nodes and the cell connection relationships are output.

[0028] S03. Numerical calculation of the flow field is performed based on the hybrid RANS-LES turbulence simulation method. Large eddy simulation is used in the flow separation prediction region, and Reynolds average method is used in the attached flow region. The output is surface pressure distribution data and velocity distribution data.

[0029] S04. The control equations are solved by coupling the WENO limiter with a high-precision shock wave capture format. The mesh is adaptively refined in the region of shock wave and boundary layer interference. The mesh node coordinates and element connection relationships of the refined mesh, as well as the updated surface pressure distribution data and velocity distribution data are output.

[0030] S05. Extract the pressure distribution data and velocity distribution data of the profile as aerodynamic performance evaluation parameters, and use the profile asymptotic optimization algorithm based on topological potential field to adjust the coordinates of the profile control points, and output the updated profile control point coordinate set.

[0031] S06. Determine whether the aerodynamic performance evaluation parameters meet the convergence criteria. If they do, output the optimized wind tunnel internal surface geometry model. If they do not meet the criteria, use the updated surface control point coordinate set as the new initial surface geometry model and return to step S02 to regenerate the mesh and calculate the flow field until convergence.

[0032] The identification criterion for high-gradient feature points is that the difference in surface curvature between adjacent control points is greater than 0.15. The surface feature curvature distribution data includes the curvature values ​​and curvature direction vectors of each control point position on the surface. The high-gradient feature point coordinate set is used to guide the determination of the mesh refinement region in the graph theory-driven mesh adaptive topology reconstruction algorithm.

[0033] The implementation of the graph theory-driven adaptive topology reconstruction algorithm is as follows: The computational grid is abstracted as a graph structure, with grid nodes corresponding to graph vertices and cell connections corresponding to graph edges. Regions with poor grid quality are identified through graph connectivity analysis. A graph cut algorithm is used to separate the sub-regions requiring reconstruction. This algorithm divides the grid domain into high-quality and low-quality regions based on the minimum-cut maximum-flow theorem. The inter-layer connections of the prism-layer grid are optimized using a shortest path algorithm to ensure the shortest connection paths between adjacent grid nodes and to avoid grid intersections. A minimum spanning tree algorithm is used to determine the priority order of grid refinement. This algorithm constructs a spanning tree using grid cell quality indices as edge weights, prioritizing the refinement of cells with the highest weights. A spectral clustering method is introduced to aggregate cells with similar flow characteristics. This method identifies regions with similar flow field structures by calculating the eigenvectors of the grid's Laplacian matrix, guiding the generation direction of anisotropic grids.

[0034] The graph theory-driven adaptive topology reconstruction algorithm employs a three-step precision approach for progressive optimization. The first step corresponds to the coarse mesh stage, where the graph structure is simplified to the backbone topology, retaining only key connections. The ratio of the mesh cell size to the wind tunnel feature length is... ∈[0.05, 0.10], quickly identifying large-scale flow separation regions. The second-order accuracy corresponds to the medium-grid stage, where the graph structure is refined to a complete topology containing all connectivity relationships, and the ratio of grid cell size to wind tunnel feature length is... ∈ [0.01, 0.05), accurately capturing shock wave location and boundary layer thickness. The third-order precision corresponds to the fine mesh stage, with the graph structure dynamically adjusted to an adaptive topology based on the local refinement according to the flow field gradient. The ratio of the mesh cell size to the wind tunnel characteristic length is... ∈[0.001, 0.01), distinguishing turbulent fluctuations and small-scale vortex structures.

[0035] Regions with poor mesh quality are identified by element distortion and element aspect ratio. Regions with poor mesh quality, defined as those with an element distortion greater than 0.7 or an element aspect ratio greater than 100, require reconstruction. The anisotropic mesh is generated along the direction of the principal strain rate of the flow field, which is calculated using the eigenvectors of the velocity gradient tensor. The initial height of the prism layer mesh is calculated based on the wall Reynolds number. The value is controlled within 1. The wall Reynolds number is calculated based on the wall friction velocity and the boundary layer thickness. The wall friction velocity is calculated using the wall shear stress and fluid density. The number of prism mesh layers depends on the boundary layer thickness and the target... The boundary layer thickness is determined using the turbulent boundary layer theory formula, and the number of prism layer mesh layers ranges from 15 to 30. The polyhedral mesh has fewer elements and better numerical accuracy than the tetrahedral mesh, with an average face number of 12 to 14 per polyhedral element. The hybrid mesh node coordinates include the three-dimensional spatial coordinates of the nodes, and the element connection relationships include the node number sequence constituting each mesh element.

[0036] The implementation method of the hybrid RANS-LES turbulence simulation method is as follows: The time-averaged flow field in the attached flow region is solved using the Reynolds-averaged Navier-Stokes equations, and the Reynolds stress term in the Reynolds-averaged Navier-Stokes equations is obtained through SST. The turbulence model is closed. Large eddy simulation (LES) is used to directly solve the large-scale turbulent structure in the flow separation prediction region, while small-scale turbulence is modeled using a subgrid model. The subgrid stress is calculated using the Smagorinsky model. The interface of the hybrid RANS-LES turbulence simulation method is determined by a mixing function, which calculates a smooth transition region based on the wall distance and the turbulence length scale.

[0037] The flow separation prediction region is identified by calculating the wall shear stress and pressure gradient parameters. A flow separation prediction region is defined as one where the wall shear stress is close to zero and the pressure gradient parameter is greater than the critical pressure gradient parameter threshold. The pressure gradient parameter is defined as follows: the pressure gradient parameter equals the partial derivative of the local pressure with respect to the flow direction, divided by the reference pressure, and multiplied by the characteristic length of the wind tunnel. The reference pressure is the incoming static pressure, and the characteristic length of the wind tunnel is the diameter of the wind tunnel test section. The critical pressure gradient parameter threshold is set to 0.05.

[0038] The hybrid RANS-LES turbulence simulation method optimizes computational efficiency using a two-step accuracy approach. In the first-step accuracy stage, the Reynolds-averaged method is used globally for preliminary solutions, identifying potential separation regions and establishing the initial flow field distribution. The ratio of the time step to the flow characteristic time is... ∈[10, 50], where the flow characteristic time is the wind tunnel characteristic length divided by the incoming velocity. In the second-order accuracy stage, the large eddy simulation method is switched to the identified flow separation prediction region, while the Reynolds average method is maintained in the remaining regions. The ratio of the time step to the flow characteristic time is... ∈[0.1, 1.0].

[0039] The SST Turbulence model combined The model's accuracy in the near-wall region and The model's stability in the far field is achieved through a mixture function that smoothly transitions between the two models. The model coefficients in the Smagorinsky model range from 0.1 to 0.2 and are dynamically adjusted based on the local grid scale and the flow Reynolds number. The profile pressure distribution data includes pressure coefficient values ​​at each control point on the profile. The pressure coefficient is calculated by dividing the difference between the local pressure and the incoming static pressure by the incoming flow pressure, where the incoming flow pressure is equal to half the incoming flow density multiplied by the square of the incoming flow velocity. The velocity distribution data includes velocity vector values ​​at each location along the test section axis.

[0040] The implementation method of the high-precision shock wave capture format coupled with the WENO limiter is as follows: The AUSM+ format is used to calculate the convective flux at the cell interface. The AUSM+ format decomposes the convective flux into convection and pressure terms, which are processed separately to avoid numerical dissipation. A fifth-order WENO limiter is used to reconstruct the flow variables at the cell interface. The WENO limiter suppresses numerical oscillations near the shock wave by weighted combination of reconstructed polynomials from multiple templates.

[0041] The criterion for adaptive mesh refinement is based on a density gradient sensor, which is defined as follows: the density gradient sensor value equals the local density gradient magnitude divided by the reference density gradient, where the reference density gradient is the incoming flow density divided by the wind tunnel characteristic length. Adaptive mesh refinement is triggered when the density gradient sensor value is greater than 0.01, and the scale of the refined mesh cells is 0.5 times that of the original mesh cells. The shock wave and boundary layer interference region is the area where the density gradient sensor value is greater than 0.01, and this region includes the compression and expansion regions before and after the shock wave.

[0042] The high-precision shock wave capture scheme coupled with the WENO limiter employs a three-step accuracy approach to improve shock wave resolution. The first-step accuracy stage uses a second-order upwind scheme for coarse calculations, quickly obtaining the approximate location of the shock wave and the basic flow field structure. The calculated residuals converge to the initial residuals. The second-order accuracy stage switches to an AUSM+ scheme coupled with a third-order WENO limiter in the shock region, improving shock capture accuracy while controlling computational cost, and the computational residual converges to a fraction of the initial residual. The third-order accuracy stage employs a fifth-order WENO limiter and adaptive mesh refinement in the shock wave and boundary layer interference region to accurately distinguish shock wave-induced boundary layer separation and reattachment phenomena, and calculates residuals that converge to the initial residuals. times.

[0043] The pressure term in the AUSM+ format is corrected using a pressure diffusion term to prevent numerical instability in the low Mach number region. The template selection for the WENO limiter is determined using a smoothness indicator, which evaluates the smoothness of the template by calculating the weighted sum of squares of the derivatives of each order of the reconstructed polynomial. The shock wave region is defined as the area where the density gradient sensor reading is greater than 0.005 and less than or equal to 0.01. The updated surface pressure and velocity distribution data have higher accuracy than those output in step S03, and are used for calculating the driving force of the surface asymptotic optimization algorithm based on the topological potential field.

[0044] The implementation of the topological potential field-based asymptotic optimization algorithm for the flow profile is as follows: A field theory model is constructed with pressure distribution as the potential energy. Each control point on the flow profile is considered a particle in the field, driven by the potential energy field generated by the local pressure gradient. The potential energy field driving force is calculated based on the local pressure gradient to drive the control point's movement. The direction of this driving force is perpendicular to the flow profile, pointing towards the direction of pressure reduction. The magnitude of this driving force is proportional to the magnitude of the pressure gradient. A topology-aware mechanism is introduced to identify flow separation and reattachment regions. This mechanism identifies the critical flow topology by analyzing the singularities and branching structures of wall friction lines. Virtual repulsive forces are added in the flow separation and reattachment regions to prevent excessive deformation of the flow profile that could lead to flow deterioration. The magnitude of these virtual repulsive forces is proportional to the second derivative of the flow profile curvature.

[0045] The surface progressive optimization algorithm based on topological potential field dynamically adjusts the potential field intensity coefficient and the radius of influence of the potential field during the iterative process. In the early stage, the potential field intensity coefficient is large and the radius of influence of the potential field is wide, realizing large-scale adjustment of the overall shape. In the later stage, the potential field intensity coefficient decreases and the radius of influence of the potential field shrinks, realizing fine optimization of local features.

[0046] The topological potential-based surface progressive optimization algorithm employs a four-step accuracy approach to achieve progressive optimization from coarse to fine. The first-step accuracy stage involves the potential field intensity coefficient... The ratio of the radius of the potential field to the total length of the profile is 1.0. The ratio of the control point movement step size to the feature size of the profile is 0.20. With a value of 0.05, the overall profile of the rapid adjustment eliminates large-scale flow separation. The potential field intensity coefficient is used in the second-step accuracy stage. The ratio of the radius of the potential field to the total length of the profile is 0.5. The ratio of the control point movement step size to the feature size of the profile is 0.10. The value is 0.02, which optimizes the pressure distribution uniformity by improving the medium-scale features of the profile. The potential field intensity coefficient is set to 0.02 in the third-step accuracy stage. The ratio of the radius of the potential field to the total length of the profile is 0.2. The ratio of the control point movement step size to the feature size of the profile is 0.05. A value of 0.005 is used to finely adjust the local curvature of the profile and eliminate small-scale flow separation bubbles. The potential field intensity coefficient is set to 0.005 for the fourth-order precision stage. The ratio of the radius of the potential field to the total length of the profile is 0.05. The ratio of the control point movement step size to the feature size of the profile is 0.02. With a value of 0.001, fine-tuning the surface details achieves optimal aerodynamic performance.

[0047] The potential energy field driving force is calculated as follows: the potential energy field driving force vector equals the negative pressure gradient vector multiplied by the potential field intensity coefficient, and then multiplied by the Gaussian weighting function of the potential field's radius of action. This Gaussian weighting function decreases exponentially with increasing distance. The second derivative of the surface curvature is calculated by taking the second derivative of the surface parameter equation with respect to the arc length parameter, reflecting the rate of change of the surface shape. The updated set of surface control point coordinates is calculated by adding the displacement of the potential energy field driving force over time to the coordinates of the prototype surface control points. The time integration step size is adaptively adjusted according to stability conditions. The surface characteristic dimension is the geometric mean of the wind tunnel contraction section inlet diameter and the test section diameter. The total surface length is the sum of the axial lengths of the wind tunnel contraction section and the diffusion section.

[0048] The aerodynamic performance evaluation parameters include a pressure distribution uniformity index and a velocity distribution stability index. The pressure distribution uniformity index is obtained by calculating the ratio of the standard deviation to the average value of the profile pressure coefficient, which is calculated based on profile pressure distribution data. The velocity distribution stability index is obtained by calculating the ratio of the maximum deviation to the average velocity of the test section's axis, which is calculated based on velocity distribution data.

[0049] The convergence criterion is that the rate of change of aerodynamic performance evaluation parameters is less than 0.5% for three consecutive iterations, and the maximum adjustment of the control point coordinates is less than 0.05% of the feature size of the aerodynamic performance. The rate of change of aerodynamic performance evaluation parameters is the difference between the aerodynamic performance evaluation parameters of the current iteration and the previous iteration, divided by the absolute value of the aerodynamic performance evaluation parameters of the previous iteration. The maximum adjustment of the control point coordinates is the maximum difference between the coordinates of corresponding control points in the updated control point coordinate set and the coordinate set before the update.

[0050] Optionally, the present invention also provides a method for implementing a wind tunnel aerodynamic optimization design system based on CFD numerical simulation using a computer. The computer is equipped with a readable storage medium that stores program instructions, which are used to execute the above-described method when the computer is run.

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

[0052] The specific implementation of step S01 involves first importing the geometric data of the wind tunnel profile and constructing a three-dimensional digital model. The profile is then represented as a spline curve using parametric curve fitting technology. Control points are then discretized along the arc length direction of the profile, with the spacing between control points adaptively adjusted according to the degree of curvature variation. More densely packed control points are set in areas of greater curvature to accurately capture the geometric features of the profile. A local curvature value is calculated for each control point. Curvature calculation is based on the first and second derivatives of the profile parametric equation with respect to the arc length parameter. The curvature gradient distribution characteristics of the profile are identified by analyzing the curvature difference between adjacent control points. When the curvature difference between adjacent control points exceeds 0.15... The region is marked as a high gradient feature region and the coordinates of the corresponding control points are recorded. These high gradient feature points represent the locations of drastic changes in the geometry of the surface, providing geometrically sensitive region information for subsequent mesh generation and flow field calculation. The final output is a dataset containing the three-dimensional coordinates of all high gradient feature points, as well as the curvature values ​​and curvature direction vector information of each control point on the complete surface.

[0053] The specific implementation of step S02 involves transforming the mesh generation problem of the computational domain into a graph theory optimization problem. First, mesh nodes are mapped to vertices of the graph, and the connectivity of mesh cells is mapped to edges, constructing an initial graph structure to represent the topological connectivity characteristics of the computational domain. A graph cut algorithm, based on the minimum cut maximum flow theorem, divides the computational domain into low-quality regions requiring reconstruction and high-quality regions to be preserved. The graph cut process uses the distortion degree and aspect ratio of computational cells as quality evaluation indicators. When the cell distortion degree is greater than 0.7 or the aspect ratio is greater than 100, it is identified as a low-quality region requiring topological reconstruction. Prismatic layer meshes are generated near the wall to accurately capture boundary layer flow characteristics. The inter-layer connectivity of the prismatic layer meshes is optimized using a shortest path algorithm. This algorithm finds the path with the minimum connection distance between adjacent layer nodes, ensuring smooth transitions between mesh layers and avoiding mesh cross-distortion. The height of the first prismatic layer is determined based on the wall Reynolds number to achieve a dimensionless wall distance. The value is controlled within 1, and the total number of prism layers is set to 15 to 30 layers to cover the boundary layer thickness. Polyhedral meshes are used in the core region instead of traditional tetrahedral meshes. Polyhedral elements have an average of 12 to 14 faces, which reduces the total number of elements and improves numerical stability at the same accuracy compared to tetrahedral meshes. The priority order of mesh refinement is determined by the minimum spanning tree algorithm. This algorithm constructs a spanning tree using the mesh element quality index as the edge weight, prioritizing the refinement of elements with the largest weight (i.e., the worst quality). Spectral clustering is introduced to aggregate mesh elements with similar flow characteristics into clusters. Spectral clustering identifies regions with similar flow field structures by calculating the eigenvectors of the Laplacian matrix formed by the mesh connectivity relationships and guides the stretching generation of anisotropic meshes along the principal strain rate direction of the flow field. The principal strain rate direction is obtained through eigenvalue decomposition of the velocity gradient tensor. The entire mesh generation process adopts a three-step accuracy progressive optimization strategy. The first step is the ratio of the mesh element size to the wind tunnel characteristic length in the coarse mesh stage. Rapidly identify large-scale flow structures between 0.05 and 0.10; the second step represents the ratio of medium-sized grid stages. Accurately captures shock waves and boundary layers between 0.01 and 0.05; the third step represents the ratio of the finer mesh stage. The turbulence detail structure is resolved between 0.001 and 0.01, and the final output is complete mesh data containing the three-dimensional coordinates of all mesh nodes and the cell connection relationships.

[0054] The specific implementation of step S03 is to use a turbulence calculation method that combines the Reynolds-averaged Navier-Stokes equations with large eddy simulation. This method employs different turbulence simulation strategies in different regions of the computational domain based on the flow characteristics to balance computational accuracy and efficiency. First, the flow field is initialized using a two-step accuracy strategy. The first step uses the Reynolds-averaged method for steady-state solution across the entire domain. The ratio of the time step to the flow characteristic time is... Between 10 and 50, a basic flow field distribution is quickly established, and potential flow separation regions are identified by analyzing wall shear stress and pressure gradient. A flow separation prediction region is defined as one where the wall shear stress is close to zero and the pressure gradient parameter exceeds 0.05. The pressure gradient parameter is defined as the ratio of the rate of change of local pressure along the flow direction normalized to the reference pressure to the characteristic length of the wind tunnel. The second step switches to large eddy simulation (LES) to directly solve the large-scale turbulent structure in the identified flow separation regions, with the time step being the ratio of the flow characteristic time. To capture transient eddy structure evolution, the model size was reduced to 0.1 to 1.0. Small-scale turbulence was modeled using the Smagorinsky subgrid model, with model coefficients dynamically adjusted within the range of 0.1 to 0.2 based on the local grid scale and Reynolds number. The attached flow region was solved using the Reynolds-averaged method to obtain the time-averaged flow field, and the Reynolds stress term was obtained using the SST method. The turbulence model is closed, and this model combines... The model's accuracy in the near-wall region and The model's stability in the far field is achieved through a mixing function that smoothly transitions between the two models. The interface between the mixed Reynolds average and large eddy simulation is determined using a mixing function based on wall distance and turbulence length scale, ensuring a smooth transition of physical quantities between different turbulence simulation methods and avoiding numerical discontinuities. After the flow field calculation converges, the pressure coefficients at each control point on the profile are extracted. The pressure coefficient is defined as the difference between the local pressure and the incoming static pressure divided by the incoming kinetic pressure, where the kinetic pressure is equal to half the product of the incoming density and the square of the incoming velocity. Simultaneously, the velocity vector distribution at each position along the test section axis is extracted, outputting complete profile pressure and velocity distribution data for aerodynamic performance evaluation.

[0055] The specific implementation of step S04 involves using a high-precision numerical scheme to accurately capture the shock wave structure in the flow field and improving computational resolution through adaptive mesh refinement. First, a three-step accuracy progressive strategy is employed to gradually improve the shock wave resolution. The first step uses a second-order upwind scheme for preliminary calculations, quickly obtaining the approximate location of the shock wave and the basic flow field structure. The computational residuals converge to the initial residuals. Once the value is multiplied by 1, the next stage can begin. The second step switches to the AUSM+ scheme in the shock region. This scheme decomposes the convection flux into convection and pressure terms, which are processed separately. The pressure term is corrected using a pressure diffusion term to prevent numerical instability in the low Mach number region. Simultaneously, a third-order WENO limiter is coupled to reconstruct the flow variables at the cell interface. The WENO limiter effectively suppresses numerical oscillations near the shock by weightedly combining the reconstructed polynomials of multiple templates. Template selection is based on a smoothness indicator, which evaluates the smoothness of the template by calculating the weighted sum of the squares of the derivatives of each order of the reconstructed polynomial. The residuals are calculated to converge to 1 / 3 of the initial residuals. The third step employs a fifth-order WENO limiter in the shock wave and boundary layer interference region and triggers adaptive mesh refinement. The mesh refinement criterion is based on a density gradient sensor, defined as the ratio of the local density gradient magnitude to the reference density gradient, where the reference density gradient is the inflow density divided by the wind tunnel characteristic length. When the density gradient sensor value is greater than 0.01, adaptive refinement is triggered, reducing the mesh cell size to 0.5 times the original size. This accurately distinguishes shock wave-induced boundary layer separation and reattachment phenomena, and the calculated residual converges to the initial residual. The refined mesh has higher spatial resolution in the shock wave and boundary layer interference region, and can capture the compression and expansion wave systems before and after the shock wave, as well as the velocity gradient distribution within the boundary layer. It outputs the node coordinates and element connection relationships of the refined mesh, as well as higher-precision surface pressure distribution data and velocity distribution data recalculated based on the new mesh.

[0056] The specific implementation of step S05 involves transforming the aerodynamic optimization problem of the profile into a field-theory-driven geometric adjustment process, constructing a topological potential field model with the profile pressure distribution as the potential energy. First, each control point on the profile is considered a particle in the potential field. Each particle is driven by the potential energy field generated by the local pressure gradient. The direction of the driving force is perpendicular to the profile and points towards the direction of pressure reduction, and the magnitude of the driving force is proportional to the pressure gradient modulus. A topological sensing mechanism is introduced to identify flow separation and reattachment regions by analyzing the singularities and branching structures of the wall friction lines. Virtual repulsive forces are added to these critical flow topological regions to prevent excessive profile deformation that could lead to flow deterioration. The magnitude of the virtual repulsive force is proportional to the second derivative of the profile curvature, which is obtained by taking the second derivative of the arc length parameter from the profile parameter equation, reflecting the rate of change of the profile shape. The potential energy field driving force is calculated by multiplying the negative pressure gradient vector by the potential field intensity coefficient and then by a Gaussian weighting function. The Gaussian weighting function decays exponentially with increasing distance, achieving spatial localization of the potential field's influence. The optimization process employs a four-step accuracy asymptotic strategy, with the first step being the potential field strength coefficient. The ratio of the radius of action of the potential field to the total length of the profile is 1.0. The ratio of the control point movement step size to the surface feature size is 0.20. To quickly adjust the overall profile of the surface to eliminate large-scale flow separation, the second step coefficient is used. The ratio is 0.5. The step size ratio is 0.10. To optimize the medium-scale characteristics to 0.02 and improve the uniformity of pressure distribution, the third-order coefficient is used. The ratio is 0.2. The step size ratio is 0.05. To finely adjust the local curvature to 0.005 and eliminate small-scale flow separation bubbles, the fourth step coefficient is used. The ratio is 0.05. The step size ratio is 0.02. To achieve optimal aerodynamic performance, fine-tuning of detailed features to 0.001 is performed. Control point positions are updated by adding the displacement of the potential energy field driving force over time to the original coordinates. The time integration step size is adaptively adjusted according to numerical stability conditions, and the updated set of control point coordinates is output for the next iteration.

[0057] The specific implementation of step S06 involves evaluating the convergence status of the current optimization iteration and determining whether the termination condition has been met. First, aerodynamic performance evaluation parameters are extracted from the updated profile pressure distribution and velocity distribution data, including pressure distribution uniformity and velocity distribution stability indices. The pressure distribution uniformity index is obtained by calculating the ratio of the standard deviation to the average value of the profile pressure coefficient, and the velocity distribution stability index is obtained by calculating the ratio of the maximum deviation of the axial velocity of the test section to the average velocity. The rate of change of the aerodynamic performance evaluation parameters between the current iteration and the previous iteration is calculated. The rate of change is defined as the absolute value of the difference between the current and previous parameter values ​​divided by the previous parameter value. Simultaneously, the maximum adjustment of the profile control point coordinates is calculated; this value is the maximum difference between the updated and unupdated control point coordinates. The convergence criterion requires that the rate of change of the aerodynamic performance evaluation parameters for three consecutive iterations be less than 0.5%, and that the maximum adjustment of the profile control point coordinates be less than 0.05% of the profile characteristic dimension. The profile characteristic dimension is defined as the geometric mean of the wind tunnel contraction section inlet diameter and the test section diameter. If both convergence criteria are met simultaneously, the optimization process is considered converged, and the current optimized wind tunnel internal surface geometry model is output as the final design result. If either criterion is not met, the updated surface control point coordinate set is used as the new initial geometry model, and the process returns to step S02 to regenerate the mesh and calculate the flow field. This process is iterated until the convergence condition is met.

[0058] It should be noted that one of the key technical ideas of this invention is a graph theory-driven adaptive topology reconstruction algorithm. This algorithm transforms the mesh generation problem into a graph connectivity optimization problem. It identifies low-quality regions requiring reconstruction using a graph cut algorithm, optimizes prism layer connections using a shortest path algorithm, determines refinement priorities using a minimum spanning tree algorithm, and guides anisotropic mesh generation using spectral clustering. Compared to traditional geometry-based mesh generation methods, the graph theory-driven method can globally optimize the mesh topology, automatically identify regions with large flow gradients for local refinement, significantly reduce the total number of mesh cells while maintaining mesh quality, improve computational efficiency, and enhance numerical stability. It exhibits particularly better adaptability when dealing with complex surfaces and regions with strong flow gradients.

[0059] The second key technical approach is a hybrid Reynolds-averaged (RA) and Large Eddy Simulation (LES) method for turbulence calculation. This method employs the computationally efficient RAS method in the attached flow region and the more accurate LES method in the flow separation prediction region. A mixing function is used to achieve a smooth transition between the different turbulence simulation methods. Compared to using the RAS method alone, the hybrid method can accurately capture the transient eddy structure evolution in the flow separation region, significantly improving the prediction accuracy of separated flows. Compared to global LES, the hybrid method significantly reduces computational costs, making it possible to complete high-precision flow field calculations within a reasonable timeframe for engineering-scale wind tunnel surface optimization, achieving an optimal balance between computational accuracy and efficiency.

[0060] The third key technical approach is a topological potential field-based asymptotic optimization algorithm for flow profiles. This algorithm transforms the flow profile optimization problem into a potential field-driven dynamic process. It drives the movement of control points by constructing a field theory model with pressure distribution as the potential energy, and introduces a topological sensing mechanism to identify critical flow structures and apply virtual repulsive forces to prevent excessive deformation. Compared to traditional gradient-based optimization methods, the topological potential field algorithm can simultaneously consider global pressure distribution and local flow topological characteristics, avoiding getting trapped in local optima. It optimizes the flow profile from coarse to fine through a four-step precision asymptotic strategy, ensuring the stability of the optimization process and improving the aerodynamic performance of the final design. It is particularly suitable for flow profile optimization in wind tunnels with complex flow separation and reattachment phenomena.

[0061] The synergistic effect of the three key technological approaches described above achieves full automation and high precision in the aerodynamic optimization of wind tunnel profiles. Graph theory-driven adaptive topology reconstruction provides a high-quality computational grid for hybrid turbulence simulation, ensuring the numerical accuracy and stability of the flow field calculations. The hybrid turbulence simulation method accurately captures complex flow phenomena near the profile, providing reliable pressure and velocity distribution data for the topology potential field optimization algorithm. The topology potential field optimization algorithm intelligently adjusts the profile geometry based on the flow field calculation results, and the optimized profile drives grid reconstruction and flow field recalculation, forming a closed-loop iteration. The three technologies work together to form a complete technical chain from grid generation to flow field solution to geometry optimization. Compared with the optimization process of traditional methods that rely on experience and manual intervention, the synergistic technology system of this invention achieves a higher level of automation and optimization efficiency, enabling wind tunnel profile designs with significantly improved aerodynamic performance at a reasonable computational cost, providing strong technical support for the development of high-performance wind tunnels.

[0062] It should be noted that this invention also solves the following technical problem: the inaccurate capture of flow details due to insufficient grid resolution in the shock wave and boundary layer interference region. Traditional grid generation methods employ uniform or geometrically proportional grid distribution strategies. In regions where the shock wave and boundary layer interact strongly, the grid scale cannot adapt to the drastic changes in the flow field gradient, resulting in significant deviations in shock wave location prediction and inaccurate capture of boundary layer separation and reattachment phenomena. This invention introduces a density gradient sensor to monitor the flow field gradient distribution in real time. When the sensor value exceeds a set threshold, it automatically triggers adaptive grid refinement, reducing the grid scale in the shock wave region to 0.5 times the original scale. Furthermore, a fifth-order WENO limiter is used to reconstruct the flow variables at the cell interface, effectively suppressing numerical oscillations near the shock wave, significantly improving shock wave capture accuracy, and accurately distinguishing boundary layer interference effects. In addition, a three-tiered accuracy progression strategy is used to quickly locate the approximate shock wave position in the coarse calculation stage and concentrate computational resources on key regions in the fine calculation stage, achieving an optimal balance between computational accuracy and efficiency, thereby solving the technical problem of inaccurate capture of flow details in the shock wave and boundary layer interference region.

[0063] It should be noted that this invention also solves the following technical problem: the inaccurate aerodynamic performance evaluation caused by insufficient accuracy of turbulence simulation in the flow separation region. Traditional computational fluid dynamics methods typically employ a uniform turbulence model across the entire flow field. While the Reynolds-averaged method (RAD) is computationally efficient, it cannot accurately predict flow separation and reattachment processes. Large eddy simulation (LES), while highly accurate, has excessive computational costs, making it difficult to apply to engineering optimization. This invention establishes a hybrid RANS-LES turbulence simulation framework. In the attached flow region, the SST turbulence model is used to solve the time-averaged flow field. In the flow separation prediction region, LES is switched to directly solve the large-scale turbulent structure. A smooth transition is achieved by determining the interface based on the wall distance and turbulence length scale using a mixing function. A dual-step accuracy calculation strategy is adopted: initially, the RAD is used globally to quickly identify potential separation regions; later, LES is used locally within the separation region. This concentrates computational resources on key areas, significantly reducing computational costs while ensuring the accuracy of flow separation prediction, thus solving the technical problem of inaccurate aerodynamic performance evaluation caused by insufficient accuracy of turbulence simulation in the flow separation region.

[0064] Specifically, the principle of this invention is as follows: The fundamental reason why this invention can solve the above-mentioned technical problems lies in establishing a two-way mapping relationship between grid topology, flow field characteristics, and profile geometry. By using graph theory, the grid structure is abstracted into a dynamically reconstructable topological graph, enabling the grid generation process to adaptively adjust according to flow field gradient information, thus avoiding the limitations of traditional fixed grids that cannot adapt to changes in the flow field. The hybrid turbulence simulation method dynamically switches the calculation model based on the flow separation prediction results, reducing the overall computational cost while ensuring the calculation accuracy in key areas, overcoming the contradiction between accuracy and efficiency in a single turbulence model. The profile optimization algorithm based on the topological potential field transforms the pressure distribution into a driving force field, causing the control point to automatically move along the direction of improving aerodynamic performance. Through a topology-aware mechanism, it identifies the critical flow state and applies constraints, avoiding flow deterioration caused by excessive profile deformation. This achieves the systematization and automation of profile adjustment, thus forming a tightly coupled closed-loop optimization system among the three components, ensuring that the optimization process converges efficiently to the global optimum.

[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 first importing the geometric data of the wind tunnel profile and constructing a three-dimensional digital model. The profile is then represented as a spline curve using parametric curve fitting technology. Control points are then discretized along the arc length direction of the profile, with the spacing between control points adaptively adjusted according to the degree of curvature variation. More densely packed control points are set in areas of greater curvature to accurately capture the geometric features of the profile. A local curvature value is calculated for each control point. The curvature calculation is based on the first and second derivatives of the profile parametric equation with respect to the arc length parameter, as specifically expressed below:

[0067] ;

[0068] In the formula, For the surface in arc length parameter The curvature value at that point, in units of ; For the surface position vector function, including Three coordinate components, in units of ; This is the arc length parameter of the profile, in units of ; Let be the first derivative of the surface position vector with respect to the arc length, and let represent the tangent direction vector of the surface, which is dimensionless. Let be the second derivative of the surface position vector with respect to the arc length, and let represent the surface curvature vector, in units of . ; This is the Euclidean norm operator used to calculate the magnitude of a vector. The surface curvature gradient distribution characteristics are identified by analyzing the curvature difference between adjacent control points. The criterion for identifying high-gradient feature points is as follows:

[0069] ;

[0070] In the formula, For the first The curvature values ​​at each control point, in units of ; For the first The curvature values ​​at each control point, in units of ; For control point serial numbers; The reference length is taken as the characteristic length of the wind tunnel, which is usually taken as the diameter of the test section, in units of... When the curvature gradient exceeds the threshold of 0.15, the region is marked as a high-gradient feature region and the coordinates of the corresponding control points are recorded. These high-gradient feature points represent the locations of drastic changes in the geometry of the surface, providing geometrically sensitive region information for subsequent mesh generation and flow field calculation. The final output is a dataset containing the three-dimensional coordinates of all high-gradient feature points, as well as the curvature values ​​and curvature direction vector information of each control point on the complete surface.

[0071] The specific implementation of step S02 involves transforming the mesh generation problem of the computational domain into a graph theory optimization problem. First, mesh nodes are mapped to vertices of the graph, and the connectivity relationships of mesh cells are mapped to edges, constructing an initial graph structure to represent the topological connectivity characteristics of the computational domain. The mesh quality evaluation index is expressed as follows:

[0072] ;

[0073] ;

[0074] In the formula, The unit twist degree is dimensionless. The maximum angle within a grid cell, in units of ; For the ideal element corresponding angle, take the following for the tetrahedral element. For hexahedral elements, take The unit is ; Pi, with a value of 3.14159; The aspect ratio of the unit is dimensionless; The longest side length of the unit, in units of ; The unit is the shortest side length of the cell. .when or Regions identified as low-quality require topology reconstruction. Prismatic layer meshes are generated near the wall to accurately capture boundary layer flow characteristics. The initial height of the prismatic layer is determined based on the wall Reynolds number, as shown below:

[0075] ;

[0076] ;

[0077] In the formula, The height of the first floor of the prism is given in units of 1. ; This is a dimensionless wall distance, with an empirical value of 1. Fluid dynamic viscosity, in units of ; Fluid density, in units of ; The wall friction speed is expressed in units of 1000 m / s. ; The wall shear stress is expressed in units of 1. Wall shear stress The wall shear rate is calculated by multiplying the wall shear rate by the fluid viscosity; the wall shear rate is the modulus of the wall normal velocity gradient. The entire mesh generation process employs a three-order accuracy progressive optimization strategy, with the mesh element size controlled as follows:

[0078] ;

[0079] In the formula, For the first The scale ratio of the grid cells in the stepped structure is dimensionless. For the first The characteristic scale of the stepped grid cells, in units of ; The reference length is taken as the characteristic length of the wind tunnel, in units of... ; Generate a step number for the grid, with a value of 1, 2, or 3; the first step. The value range is from 0.05 to 0.10, second step. The value range is from 0.01 to 0.05, third step. The value ranges from 0.001 to 0.01. The final output contains complete mesh data including the 3D coordinates of all mesh nodes and the element connection relationships.

[0080] The specific implementation of step S03 is to use a turbulence calculation method that combines the Reynolds-averaged Navier-Stokes equations with large eddy simulation. This method employs different turbulence simulation strategies in different regions of the computational domain based on flow characteristics to balance computational accuracy and efficiency. The criteria for identifying the flow separation prediction region are expressed as follows:

[0081] ;

[0082] In the formula, The pressure gradient parameter is dimensionless. This is the partial derivative of pressure with respect to the flow direction coordinate, in units of... ; The reference length is taken as the characteristic length of the wind tunnel, in units of... ; The reference pressure is the incoming static pressure, and the unit is . ; Flow direction coordinates, in units of ; Local pressure, unit: When the wall shear stress Approaching zero and The time step is determined to be within the flow separation prediction region. The time step is controlled as follows:

[0083] ;

[0084] In the formula, For the first The time step of the ladder, in units of ; For the first The time step ratio coefficient of the step is dimensionless. The reference length is taken as the characteristic length of the wind tunnel, in units of... ; The incoming flow velocity is expressed in units of... ; The step number is used for flow field calculation, with a value of 1 or 2; the first step. The value range is from 10 to 50, second step. The value ranges from 0.1 to 1.0. The pressure coefficient is defined as follows:

[0085] ;

[0086] In the formula, The pressure coefficient is dimensionless. Local pressure, unit: ; The static pressure of the incoming flow is expressed in units of 1. ; The incoming flow density is expressed in units of... ; The incoming flow velocity is expressed in units of... After the flow field calculation converges, the pressure coefficients at each control point on the surface and the velocity vector distribution at each position on the test section axis are extracted, and complete surface pressure distribution data and velocity distribution data are output for aerodynamic performance evaluation.

[0087] The specific implementation of step S04 involves using a high-precision numerical format to accurately capture the shock wave structure in the flow field and improving computational resolution through adaptive mesh refinement. The mesh refinement criterion is based on a density gradient sensor and is specifically expressed as follows:

[0088] ;

[0089] In the formula, It is a density gradient sensor, dimensionless; This is the density gradient vector, in units of ; This is the Euclidean norm operator, used to calculate the magnitude of a vector. The incoming flow density is expressed in units of... ; The reference length is taken as the characteristic length of the wind tunnel, in units of... ; Local density, in units of .when Adaptive encryption is triggered at any time, and the mesh element scale is adjusted as follows after encryption:

[0090] ;

[0091] In the formula, This refers to the scale of the refined grid cells, in units of... ; This is the original grid cell scale, in units of A three-step accuracy progressive strategy is adopted to gradually improve the shock wave resolution. The residual convergence criterion is expressed as follows:

[0092] ;

[0093] In the formula, For the first The residual ratio of the steps is dimensionless. For the first The flux residual vector of the next iteration, in units of ; This is the initial flux residual vector, in units of ; This represents the number of iterations. The shock wave capture step number is set, with a value of 1, 2, or 3; the first step converges to... The second step converges to The third step converges to The refined mesh has higher spatial resolution in the shock wave and boundary layer interference region, and can capture the compression wave system and expansion wave system before and after the shock wave, as well as the velocity gradient distribution within the boundary layer. It outputs the node coordinates and element connection relationships of the refined mesh, as well as higher-precision surface pressure distribution data and velocity distribution data recalculated based on the new mesh.

[0094] The specific implementation of step S05 involves transforming the aerodynamic optimization problem of the profile into a field-theory-driven geometric adjustment process, constructing a topological potential field model with the profile pressure distribution as the potential energy. The calculation method for the driving force of the potential energy field is as follows:

[0095] ;

[0096] In the formula, For the first The direction vector of the potential energy field driving force on each control point is dimensionless. For the first The potential field strength coefficient of the step is dimensionless. For the first Pressure coefficient gradient vector at each control point, in units of ; The magnitude of the pressure coefficient gradient vector, in units of . ; For the first The position vector of each control point, in units of ; For the first The position vectors of adjacent control points, in units of ; The adjacent control point number; For the first The set of adjacent control points of each control point; For the first The ratio of the potential field radii of the steps is dimensionless. The total length of the profile, in units of ; To optimize the ladder sequence number, the value can be 1, 2, 3 or 4; This represents the current control point number. Gaussian weighting function. The exponential term is the ratio of the squares of the two lengths, and therefore dimensionless. This function decays exponentially with increasing distance, realizing the spatial localization of the potential field. The virtual repulsive force is calculated as follows:

[0097] ;

[0098] In the formula, For the first The virtual repulsive force direction vector at each control point is dimensionless. This is the repulsion coefficient, with an empirical value of 0.1; For the first The second derivative of the surface curvature at each control point, in units of ; As the reference curvature, the value is [value to be filled in]. The unit is ; The reference length is taken as the characteristic length of the wind tunnel, in units of... ; For the first The surface normal unit vector at each control point is dimensionless. This is the arc length parameter, in units of The unit vector normal to the surface. Perpendicular to the tangential direction of the profile, it is obtained by normalizing the cross product of the profile tangential vector and the reference direction vector. The tangential vector is calculated using the first derivative of the profile position vector with respect to the arc length. The control point position update is represented as follows:

[0099] ;

[0100] In the formula, For the updated number The position vector of each control point, in units of ; For the previous version The position vector of each control point, in units of ; For the first The ratio of control point movement steps in a stepped structure, dimensionless; These are the surface feature dimensions, in units of ; The direction vector of the potential energy field driving force is dimensionless. This is a dimensionless virtual repulsive force direction vector. (Surface feature dimensions) It is expressed as follows:

[0101] ;

[0102] In the formula, These are the surface feature dimensions, in units of ; The diameter of the wind tunnel's contraction section entrance, in units of ; The diameter of the test section is given in units of 1. The optimization process employs a four-step accuracy progressive strategy, with the first step... Second step Third step Fourth step The updated set of control point coordinates is output for the next iteration.

[0103] The specific implementation of step S06 involves evaluating the convergence state of the current optimization iteration and determining whether the termination condition has been met. Aerodynamic performance evaluation parameters include pressure distribution uniformity indices and velocity distribution stability indices, specifically represented as follows:

[0104] ;

[0105] ;

[0106] In the formula, It is a dimensionless index for pressure distribution uniformity. is the standard deviation of the surface pressure coefficient, which is dimensionless; This represents the average value of the surface pressure coefficient, which is dimensionless. It is a dimensionless index for the stationarity of velocity distribution; The maximum deviation of the axis speed of the test section, in units of ; The average speed of the test segment axis, in units of Standard deviation of the surface pressure coefficient It is expressed as follows:

[0107] ;

[0108] In the formula, is the standard deviation of the surface pressure coefficient, which is dimensionless; This represents the total number of surface control points. For the first Pressure coefficient at each control point, dimensionless; This represents the average value of the surface pressure coefficient, which is dimensionless. Control point number. Average value of the profile pressure coefficient. It is expressed as follows:

[0109] ;

[0110] In the formula, This represents the average value of the surface pressure coefficient, which is dimensionless. This represents the total number of surface control points. For the first Pressure coefficient at each control point, dimensionless; Here are the control point numbers. The rate of change of aerodynamic performance evaluation parameters is expressed as follows:

[0111] ;

[0112] In the formula, For the first The rate of change of aerodynamic performance evaluation parameters in each iteration is dimensionless. For the first The aerodynamic performance evaluation parameter values ​​for each iteration are dimensionless. For the first The aerodynamic performance evaluation parameter values ​​for each iteration are dimensionless. The iteration count is given. The maximum adjustment amount of the surface control point coordinates is expressed as follows:

[0113] ;

[0114] In the formula, This represents the maximum adjustment amount of the control point coordinates for the profile, in units of... ; For the updated number The position vector of each control point, in units of ; For the previous version The position vector of each control point, in units of ; This represents the total number of surface control points. Here are the control point indices. The convergence criterion requires that the following conditions be met in three consecutive iterations. and ,in These are the surface feature dimensions, in units of If both convergence criteria are met simultaneously, the optimization process is considered converged, and the current optimized wind tunnel internal geometry model is output as the final design result. If either criterion is not met, the updated set of control point coordinates is used as the new initial geometry model, and the process returns to step S02 to regenerate the mesh and calculate the flow field. This process is iterated until the convergence condition is met.

[0115] It should be explained that the potential energy field driving force formula in step S05 is based on the gradient descent principle. It guides the adjustment of the orientation towards the low-pressure region through the negative pressure gradient direction, thereby adjusting the pressure coefficient gradient vector. Divide by its modulus The unit direction vector is obtained to achieve dimensionlessness. This unit vector is multiplied by the summation term of the Gaussian weighting function and then multiplied by the potential field intensity coefficient. Obtain the dimensionless driving force direction vector Gaussian weighting function The numerator of the exponential term is the square of the distance between control points, in units of... The denominator is the square unit of the potential field's range of action. Since both have the same dimensions, the exponent term becomes dimensionless. This function achieves spatial localization of the potential field, avoiding global disturbances, and includes the set of adjacent control points. By determining that the distance between control points is less than the radius of the potential field. To determine this. The virtual repulsion formula is derived from the second derivative of curvature. Units are Divided by reference curvature Units are We obtain a dimensionless ratio, which is proportional to the normal unit vector. and repulsion coefficient Multiplying them yields a dimensionless virtual repulsive force direction vector. To prevent excessive deformation of the flow separation region profile, virtual repulsion coefficient The empirical value is 0.1, which can be adjusted according to specific working conditions. The control point position update formula will use the dimensionless driving force vector. With dimensionless virtual repulsive vector The sums are used to obtain the total dimensionless adjustment direction vector, which is then multiplied by the feature dimension of the profile. Units are Ratio of moving step size The dimensionless actual displacement vector unit is obtained as Add the displacement vector to the original position vector Units are The new position vector is obtained. Units are The dimensions on both sides of the formula are Maintain consistency. The four-step precision asymptotic strategy gradually reduces the potential field strength coefficient. The ratio of the effective radius was reduced from 1.0 to 0.05. Decrease the step size ratio from 0.20 to 0.02. Reducing the value from 0.05 to 0.001 achieves an optimization path from the overall profile to local details. This strategy transforms the complex aerodynamic optimization problem into a field theory-driven geometric adjustment process, effectively improving the optimization convergence speed and avoiding getting trapped in local optima. This results in efficient and accurate optimization of the aerodynamic performance of the wind tunnel profile. The convergence criterion formula in step S06 uses the rate of change of aerodynamic performance evaluation parameters. Ratio of the maximum adjustment amount of the control point coordinates of the profile Two dimensionless indices are used to determine the convergence state of the optimization: the pressure distribution uniformity index. Standard deviation Compared with the average The ratios are all dimensionless, and the velocity distribution stability index For the maximum deviation Units are With average speed Units are The ratio is also dimensionless. Both indicators are dimensionless, which facilitates comparison between different working conditions. The constraint of three consecutive iterations ensures the stability and reliability of convergence. The setting of the change rate threshold of 0.005 (0.5%) and the adjustment ratio threshold of 0.0005 (0.05%) balances the calculation accuracy and efficiency. This formula effectively avoids insufficient optimization caused by premature termination and waste of computational resources caused by excessive iteration, and provides a reliable termination criterion for the aerodynamic optimization of the wind tunnel surface.

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

[0117] Table 1. Variable Explanation Table (Part 1)

[0118]

[0119] Table 2. Variable Explanation Table (Part Two)

[0120]

[0121] To better understand and implement this invention, a specific application scenario is provided below as Example 2: A technical team undertook the task of optimizing the internal profile of a transonic wind tunnel. The test section of the wind tunnel has a diameter of 1.2m, the inlet diameter of the contraction section is 3.6m, the design Mach number is 0.85, the total inflow pressure is 101325Pa, and the total inflow temperature is 288K. The original wind tunnel internal profile exhibited significant pressure fluctuations and uneven velocity distribution in the test section during high-speed operation, affecting the accuracy of the experimental data. The technical team decided to redesign the wind tunnel internal profile using the aerodynamic optimization design method of this invention.

[0122] The technical team first established an initial geometric model of the wind tunnel's profile, with a total length of 8.4m for the contraction and diffusion sections. The profile was defined by 68 control points. By extracting the characteristic curvature distribution of the profile, they discovered regions of drastic curvature changes in the latter half of the contraction section and at the entrance of the test section. The maximum curvature difference between adjacent control points reached 0.28. It far exceeds 0.15 The high gradient feature point identification criteria were established. The technical team marked 15 high gradient feature points in these areas, and their coordinate data are shown in Table 3.

[0123] Table 3. Coordinate data of high gradient feature points

[0124]

[0125] The technical team employed a graph theory-driven adaptive topology reconstruction algorithm to generate hybrid meshes. In the first-order accuracy stage, the ratio of mesh cell size to wind tunnel feature length... Setting the ratio to 0.08 generated a coarse mesh containing 1.85 million elements, quickly identifying the large-scale flow separation region in the latter half of the contraction section. In the second-order accuracy stage, the ratio... With a setting of 0.03, the number of grid cells increased to 5.62 million, accurately capturing the shock wave location at the test section entrance at an axial coordinate of 4.15m. In the third-order accuracy stage, the ratio... The mesh size was set to 0.005, and local refinement was performed near the shock wave and in the boundary layer region, resulting in a final mesh count of 12.48 million cells. A 22-layer prism mesh was generated near the wall, with the first layer having a height of 0.008 mm to ensure... The value is controlled within 0.85. The core region uses a polyhedral mesh, with an average of 13 faces per polyhedral unit, which reduces the number of units by about 40% compared to the traditional tetrahedral mesh.

[0126] The technical team performed numerical calculations of the flow field based on the hybrid RANS-LES turbulence simulation method. In the first-order accuracy stage, SST was used across the entire domain. The Reynolds-averaged solution for the turbulence model is obtained by the ratio of the time step to the characteristic flow time. The time step was set to 25, and the flow characteristic time was 0.0107 s. After 500 time steps, the region from 2.8 m to 3.5 m along the axial coordinates of the latter half of the contraction section was identified as a potential flow separation region. In this region, the wall shear stress decreased to 0.15 Pa, and the pressure gradient parameter reached 0.078, exceeding the critical threshold of 0.05. In the second-order accuracy stage, the large eddy simulation method was switched to the identified flow separation prediction region. The Smagorinsky model coefficients were set to 0.15, and the ratio of the time step to the flow characteristic time was... The value was set to 0.5, and the Reynolds average method was maintained for the remaining regions. After 2000 time steps of calculation, the surface pressure distribution data and the axial velocity distribution data of the test section were obtained, such as... Figure 2 As shown.

[0127] The technical team employed a high-precision shock wave capture scheme coupled with a WENO limiter to solve the governing equations. In the first-order accuracy stage, a second-order upwind scheme was used for coarse calculations, and the residuals were calculated from the initial values. convergence to This allows for the rapid acquisition of the approximate shock wave location. In the second-order accuracy stage, the shock wave region at the test section inlet axial coordinates from 4.05m to 4.25m is switched to the AUSM+ scheme coupled with a third-order WENO limiter, and the calculated residuals converge to... The density gradient sensor reaches its maximum value of 0.032 at the shock wave location, such as... Figure 3As shown. In the third-order accuracy stage, a fifth-order WENO limiter is used in the shock wave and boundary layer interference region, and adaptive mesh refinement is performed, reducing the mesh cell size in this region to 0.5 times the original size and adding 1.58 million new mesh cells. The calculated residuals eventually converge to... It accurately distinguishes the boundary layer separation and reattachment phenomena induced by shock waves.

[0128] The technical team employed a topological potential-based progressive optimization algorithm to adjust the coordinates of the profile control points. In the first-order accuracy stage, the potential field intensity coefficient... Set to 1.0, the ratio of the radius of the potential field to the total length of the profile. Set to 0.20, the ratio of the control point movement step size to the feature size of the profile. The potential energy coefficient was set to 0.05, and the surface feature size was 2.08m. By constructing a pressure distribution potential energy field, the driving force of the potential energy field at control point P3 was calculated to be the maximum, reaching 125N, driving the control point to move radially outward by 0.104m. After 8 iterations, the large-scale flow separation in the latter half of the contraction section was successfully eliminated, and the pressure distribution uniformity index decreased from the initial 0.182 to 0.095. In the second-step accuracy stage, the potential field strength coefficient... Set to 0.5, ratio Set to 0.10, ratio Set to 0.02, after 15 iterations, the medium-scale features of the surface were optimized, and the stability index of the axial velocity distribution in the test section decreased from 0.068 to 0.028. In the third-order accuracy stage, the potential field intensity coefficient... Set to 0.2, ratio Set to 0.05, ratio Set to 0.005, and after 23 iterations, the local curvature of the surface was finely adjusted, eliminating small-scale flow separation bubbles at the inlet of the test section. In the fourth-order accuracy stage, the potential field intensity coefficient... Set to 0.05, ratio Set to 0.02, ratio The value was set to 0.001, and after 12 iterations, the surface details were fine-tuned to achieve optimal aerodynamic performance.

[0129] The technical team introduced a topology-aware mechanism to identify flow separation and reattachment regions during the optimization process. By analyzing the singularities and branching structures of the wall friction lines, three critical flow topology regions were identified, located at axial coordinates of 2.95m, 3.48m, and 4.12m, respectively. Virtual repulsive forces were added to these regions to prevent excessive deformation of the profile. The magnitude of the virtual repulsive force is proportional to the second derivative of the profile curvature, with a maximum virtual repulsive force of 85N, effectively avoiding uneven and undesirable profile shapes.

[0130] After a total of 58 iterations, the aerodynamic performance evaluation parameters met the convergence criteria. The change rates of the pressure distribution uniformity index in the three consecutive iterations were 0.38%, 0.32%, and 0.28%, respectively, and the change rates of the velocity distribution stability index were 0.42%, 0.35%, and 0.31%, respectively, all less than 0.5%. The maximum adjustment of the profile control point coordinates decreased from 0.0012m in the 55th iteration to 0.00085m in the 58th iteration, which is less than 0.05% of the profile feature size, i.e., 0.00104m. The technical team output the optimized wind tunnel profile geometry model. The optimized profile has a pressure distribution uniformity index of 0.015 and a velocity distribution stability index of 0.012, representing reductions of 91.8% and 82.4%, respectively, compared to the initial profile.

[0131] The advancements of this invention over traditional methods are mainly reflected in three aspects. First, the graph-driven adaptive topology reconstruction algorithm abstracts the mesh into a graph structure and optimizes the mesh topology using graph connectivity analysis and shortest path algorithms, avoiding the problems of poor mesh quality and chaotic connectivity in traditional mesh generation methods. Simultaneously, the three-step accuracy method achieves progressive mesh refinement from coarse to fine, significantly improving mesh generation efficiency and quality. Second, the hybrid RANS-LES turbulence simulation method combines the efficiency of the Reynolds-averaged method with the high accuracy of the Large Eddy Simulation method. Through a mixing function, it achieves a smooth transition between the attached flow region and the separated flow region, overcoming the limitations of insufficient accuracy or excessive computational cost of traditional single turbulence models in complex flows. The two-step accuracy method further optimizes computational efficiency. Third, the topological potential field-based progressive optimization algorithm transforms the pressure distribution into a potential energy field. The potential energy field drives the deformation direction of the profile automatically. The topological sensing mechanism is introduced to identify the critical flow region and apply virtual repulsive force constraints. This avoids the subjectivity and blindness of manually setting optimization directions and constraints in traditional optimization methods. The four-step precision method realizes the step-by-step optimization from the overall profile to local details, ensuring the stability and convergence of the optimization process.

[0132] 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 method for aerodynamic optimization design of wind tunnel surfaces, characterized in that, Establish an initial geometric model of the wind tunnel surface and extract the surface curvature distribution. Mark high gradient feature points in areas of drastic curvature change. A graph-theory-driven adaptive topology reconstruction algorithm is used to generate hybrid meshes, generating prism-layer meshes near the walls and polyhedral meshes in the core region. Numerical calculations of the flow field were performed based on the hybrid RANS-LES turbulence simulation method. Large eddy simulation was used in the flow separation prediction region, while Reynolds average method was used in the attached flow region. A high-precision shock wave capture scheme coupled with a WENO limiter is used to solve the control equations, and adaptive mesh refinement is performed in the shock wave and boundary layer interference region. Pressure and velocity distribution data of the wind tunnel profile are extracted as aerodynamic performance evaluation parameters, and the coordinates of the control points of the wind tunnel profile are adjusted using a topological potential field-based asymptotic optimization algorithm. It is then determined whether the aerodynamic performance evaluation parameters meet the convergence criteria. If they do, the optimized wind tunnel profile geometry model is output; otherwise, the process is repeated until convergence.

2. The method according to claim 1, characterized in that, The identification criterion for high gradient feature points is that the difference in surface curvature between adjacent control points is greater than 0.15 rad / m.

3. The method according to claim 2, characterized in that, The graph theory-driven adaptive topology reconstruction algorithm specifically abstracts the computational grid into a graph structure, where grid nodes correspond to graph vertices and cell connection relationships correspond to graph edges, and identifies regions with poor grid quality through graph connectivity analysis.

4. The method according to claim 3, characterized in that, In the graph theory-driven adaptive topology reconstruction algorithm, the graph cut algorithm is used to separate the sub-regions that need to be reconstructed, the inter-layer connectivity of the prism layer grid is optimized based on the shortest path algorithm, the minimum spanning tree algorithm is used to determine the priority order of grid densification, and the spectral clustering method is introduced to aggregate units with similar flow characteristics.

5. The method according to claim 4, characterized in that, The graph theory-driven adaptive topology reconstruction algorithm employs a three-tiered precision approach for progressive optimization. The first tier of precision corresponds to the ratio of the grid cell size to the wind tunnel feature length in the coarse grid stage (∈ [0.05, 0.10]), the second tier of precision corresponds to the ratio of the grid cell size to the wind tunnel feature length in the medium grid stage (∈ [0.01, 0.05)), and the third tier of precision corresponds to the ratio of the grid cell size to the wind tunnel feature length in the fine grid stage (∈ [0.001, 0.01)).

6. The method according to claim 5, characterized in that, The height of the first layer of the prism grid is calculated based on the wall Reynolds number. The value is controlled within 1, and the number of prism mesh layers ranges from 15 to 30 layers.

7. The method according to claim 6, characterized in that, In the hybrid RANS-LES turbulence simulation method, the time-averaged flow field in the attached flow region is solved using the Reynolds-averaged Navier-Stokes equations. The Reynolds stress term in the Reynolds-averaged Navier-Stokes equations is obtained through SST. The turbulence model is closed, and large eddy simulation is used to directly solve the large-scale turbulence structure in the flow separation prediction region, while small-scale turbulence is modeled through a subgrid model.

8. The method according to claim 7, characterized in that, The method for identifying the flow separation prediction region is to calculate the wall shear stress and pressure gradient parameters. When the wall shear stress is close to zero and the pressure gradient parameter is greater than the critical pressure gradient parameter threshold, it is determined to be a flow separation prediction region. The critical pressure gradient parameter threshold is 0.

05.

9. The method according to claim 8, characterized in that, The hybrid RANS-LES turbulence simulation method adopts a two-step accuracy approach to optimize computational efficiency. In the first-step accuracy stage, the Reynolds average method is used for preliminary solution across the entire domain, with the ratio of the time step to the flow characteristic time ∈ [10, 50]. In the second-step accuracy stage, the method is switched to large eddy simulation in the identified flow separation prediction region, with the ratio of the time step to the flow characteristic time ∈ [0.1, 1.0].

10. The method according to claim 9, characterized in that, The high-precision shock wave capture format coupled with the WENO limiter specifically uses the AUSM+ format to calculate the convection flux of the cell interface and a fifth-order WENO limiter to reconstruct the flow variables of the cell interface.

Images