Suborbital vehicle performance optimization method considering flight state constraints

Abstract

The application discloses a sub-orbital vehicle carrying performance optimization method considering flight state constraints and belongs to the field of sub-orbital vehicle design optimization.The application realizes the method as follows: a multi-disciplinary analysis model and a multi-disciplinary design optimization problem model of a sub-orbital vehicle are established by comprehensively considering the coupling relationship among key disciplines such as geometry and aerodynamics of the sub-orbital vehicle;an approximate optimization strategy based on an initial sample expansion mechanism is adopted to optimize design variables of the sub-orbital vehicle with the maximum effective payload of the sub-orbital vehicle as an objective function, considering flight state constraints such as overload and dynamic pressure;in the optimization process, a Kriging surrogate model is adopted to replace a high-time-consuming analysis model of the sub-orbital vehicle, a constraint differential evolution algorithm is used to solve a pseudo-optimal solution of the surrogate model, and the surrogate model is managed and updated by combining a comprehensive improvement probability and a key design space method, so that the sub-orbital vehicle carrying performance optimization process is efficiently guided to converge to a global optimal solution, that is, the sub-orbital vehicle carrying performance optimization is realized.

Description

TECHNICAL FIELD

[0001] The application relates to a suborbital vehicle carrying performance optimization method considering flight state constraints and belongs to the technical field of suborbital vehicle design optimization. BACKGROUND

[0002] A suborbital vehicle is a new type of vehicle that can reach the space boundary (about 100 km). Due to its reusability, maintainability and high efficiency, the suborbital vehicle has attracted widespread attention from the tourism industry, transportation industry and other industries (Chang, 2020). In recent years, suborbital vehicles such as Spaceship One, SpaceLiner and Lynx have been fully developed. However, the suborbital vehicle often faces challenges such as strong overload, high heat flux and high dynamic pressure during its re-entry and return to the earth's atmosphere. In order to cope with the above challenges, domestic researchers have done a lot of research on the optimization of suborbital vehicle subsystem design. Wang (2018) uses the direct search region method to optimize the thermal protection system of the vehicle, so as to minimize the stagnation point heat flux during flight. Zhang (2014) designs an optimal attack angle program for the re-entry and return process of the suborbital vehicle based on the particle swarm optimization algorithm. Luo (2023) establishes a RBCC engine performance model based on the quasi-one-dimensional thermodynamic method, and explores the specific impulse increase ratio and the thrust-area ratio Pareto front based on the NSGA-II method. However, the overall performance of the suborbital vehicle is affected by multiple subsystems such as geometric characteristics and aerodynamics, so the suborbital vehicle design optimization is a typical Multidisciplinary Design Optimization (MDO) problem. In addition, the suborbital vehicle design optimization often involves high-time-consuming analysis models such as aerodynamic force and heat. If the traditional evolutionary optimization algorithm is still used to solve the suborbital vehicle design optimization problem, the computational cost will increase dramatically.

[0003] In order to improve the design efficiency of the overall scheme stage of the suborbital vehicle, it is necessary to develop a suborbital vehicle carrying performance optimization method considering flight state constraints, so as to improve the design optimization efficiency and shorten the design cycle, so as to quickly realize the rapid optimization and modification of the suborbital vehicle design scheme in the overall design stage, and provide scientific basis and reference for the suborbital vehicle system scheme demonstration and overall design.

[0004] In order to better illustrate the technical scheme of the application, the mathematical basis related to the Kriging surrogate model is introduced as follows.

[0005] The Kriging surrogate model is an unbiased optimal estimation interpolation model, which is composed of a global model and a local bias term, as shown in formula (1)

[0006]

[0007] where μ(x) is a polynomial global approximation model in the design space, and μ is usually a constant when the value of the approximated object is unknown; the local deviation term Z(x) is a stochastic process with zero mean and non-zero variance σ 2 , and its covariance is non-zero, which represents the local deviation based on the global approximation model. The correlation function and the covariance matrix can be represented by equation (2) and equation (3), respectively.

[0008]

[0009] Cov[Z(x i ),Z(x j )]=σ 2 R[R(x i ,x j )] (3)

[0010] where R is a symmetric correlation matrix; R(·,·) is a correlation function; n v is the dimension of the design variable; x i is the training sample point. The estimated values of μ and σ 2 can be obtained by the least square method according to equation (4):

[0011]

[0012] where n KRG represents the number of sample points for constructing the surrogate model; I represents an n KRG -dimensional unit row vector. The correlation coefficient θ k

[0013]

[0014] In addition, the Kriging surrogate model can estimate the variance s 2 (x) of the predicted value at any point x, so as to evaluate the approximation error of the surrogate model. s 2 (x) can be represented as:

[0015]

[0016] where :r represents a correlation function vector, and its specific form is as follows

[0017] SUMMARY

[0018] To address the challenge of improving the payload performance of suborbital vehicles, this invention aims to provide a method for optimizing suborbital vehicle payload performance considering flight state constraints. This method guides the suborbital vehicle design optimization process towards the optimal solution rapidly, while ensuring that the suborbital vehicle design meets constraints such as overload, thermal flux, and dynamic pressure, thereby increasing the vehicle's effective payload. This invention is applicable to the overall design optimization of suborbital vehicles with diverse mission requirements, improving both the performance and design efficiency of suborbital vehicles.

[0019] The objective of this invention is achieved through the following technical solution.

[0020] This invention discloses a method for optimizing the launch performance of suborbital vehicles considering flight state constraints. It comprehensively considers the coupling relationships between key disciplines such as geometry and aerodynamics of suborbital vehicles, establishing a multi-disciplinary analysis model and a multi-disciplinary design optimization problem model for suborbital vehicles. This invention employs an approximate optimization strategy based on an initial sample expansion mechanism, using the maximum effective payload of the suborbital vehicle as the objective function, and considering flight state constraints such as overload and dynamic pressure to optimize the design variables of the suborbital vehicle. During the optimization process, a Kriging surrogate model is used to replace the time-consuming analysis model of the suborbital vehicle, and the pseudo-optimal solution of the surrogate model is solved by a constrained differential evolution algorithm. Combining the comprehensive improvement probability and key design space methods, the surrogate model is managed and updated, efficiently guiding the suborbital vehicle launch performance optimization process to converge to the global optimum, thus achieving suborbital vehicle launch performance optimization considering flight state constraints. This method overcomes the problems of low solution efficiency and poor global convergence of traditional optimization methods, enabling efficient discovery of high-performance design schemes for suborbital vehicles and shortening the design cycle.

[0021] The suborbital spacecraft launch performance optimization method considering flight state constraints disclosed in this invention includes the following steps:

[0022] Step 1: Determine the suborbital spacecraft launch performance optimization model and determine the initial parameters of the approximate optimization strategy based on the initial sample expansion mechanism.

[0023] The specific implementation method for step one is as follows:

[0024] With suborbital spacecraft payload mass m pl The objective function is to maximize the takeoff mass m. takeoff Axial overload N a Normal overload N n Reentry stage end velocity V end Reentry phase range R, reentry phase end lift-to-drag ratio C l / C d Flight kinetic pressure Q, terminal velocity angle Θ during reentry touch stagnation heat flux density qstag and wing leading edge heat flux density q wing Constraints, establish suborbital vehicle suborbital vehicle carrying performance optimization model, as shown in equation (8):

[0025] find X=[c v,1 ,c v,2 ,c h,1 ,c h,2 ,L int ,L ext ,b root ,b tip ,χ int ,χ ext ,L h ,W r ,f u ,f a ,k α ]

[0026] min f(X)=-m pl

[0027]

[0028] In the formula, X represents the design variable, specifically including vehicle head profile vertical direction control point 1 coordinate c v,1 , vehicle head profile vertical direction control point 2 coordinate c v,2 , vehicle head profile horizontal direction control point 1 coordinate c h,1 , vehicle head profile vertical direction control point 2 coordinate c h,2 , inner wing segment half span length L int , outer wing segment half span length L ext , wing root chord length b root , wing tip chord length b tip , inner wing segment sweepback angle χ int , outer wing segment sweepback angle χ ext , head length L h , afterbody width W r , unpredictable fuel mass ratio f u , additional fuel mass ratio f a and angle of attack reduction rate k α as design variables. X LB and X UB respectively represent the lower and upper limits of the design variable.

[0029] The initial parameters of the approximate optimization strategy based on the initial sample expansion mechanism include the initial sampling number n doe , the constraint threshold scaling factor η, the initial sample expansion number n l , the number of newly added sample points in the local search stage n a and the maximum model call number

[0030] Step two: based on the suborbital vehicle design variables, a parameterized geometric master model of the suborbital vehicle is established, and the geometric characteristics of the suborbital vehicle are obtained. The geometric characteristics of the suborbital vehicle include the reference area, the maximum curvature radius of the head point, the minimum curvature radius of the head point, the curvature radius of the wing leading edge, the sweep angle of the wing leading edge, the total volume, the total surface area, the wing area and the tail area.

[0031] The specific implementation method of step two is as follows:

[0032] Based on the NURBS method shown in formula (9), the head, the fuselage, the wing, the afterbody and the tail of the suborbital vehicle are parameterized modeled to obtain the parameterized geometric master model:

[0033]

[0034] In the formula, d i is the curve control point, R i,p (u) is the rational basis function defined on u∈[0, 1], ω i is the weight coefficient, N i,p (u) is the i-th p-order B-spline basis function, and n represents the number of nodes. N i,p (u) can be represented by formula (10) and formula (11); when the node u and the order p are determined, the B-spline basis function can be uniquely determined

[0035]

[0036] In the formula, N i,p (u) is the i-th p-order B-spline basis function, and u i is the i-th node.

[0037] The parameterized shape of the head, the fuselage, the wing, the afterbody and the tail is generated by formula (9)-(11); and the geometric characteristics (reference area, maximum curvature radius of head point, minimum curvature radius of head point, curvature radius of wing leading edge, sweep angle of wing leading edge, total volume, total surface area, wing area and tail area) are measured as inputs for obtaining the lift coefficient, the drag coefficient, the head point heat flux density, the wing inner wing segment leading edge heat flux density, the wing outer wing segment leading edge heat flux density, the reentry return mass and the trajectory model.

[0038] Step three: according to the geometric master model of the suborbital vehicle and the reference area in step two, the lift coefficient and the drag coefficient of the suborbital vehicle are obtained.

[0039] The specific implementation method of step three is as follows:

[0040] ①. Based on the parameterized geometry main model and the computational domain in step two, unstructured mesh is divided; and local refinement is added at the head stagnation point and the wing.

[0041] Boundary layer mesh is added at the head and the wing, and the first layer mesh height is shown in equation (12)

[0042]

[0043] In the equation, dimensionless wall distance y + ≈1、L r is the reference length of the aircraft, and Re is the Reynolds number, which is represented by equation (13)

[0044]

[0045] In the equation, p ∞ is the atmospheric density, v ∞ is the incoming flow speed in the flight condition, and m is the atmospheric viscosity.

[0046] After local refinement and the addition of boundary layer mesh, mesh A is obtained.

[0047] ②. Based on the reference area of the sub-orbital aircraft measured in step two and mesh A in ①, the lift coefficient and the drag coefficient of the sub-orbital aircraft are solved by the CFD method. In the CFD solution, the density-based solver and the energy equation are used, and the incoming flow is set as the pressure far-field boundary condition; the spatial discretization is performed using the upwind splitting format, and the Green-Gauss node method is used to calculate the gradient; in addition, 3-level FMG is used for initialization before starting the solution to improve the solution convergence speed.

[0048] Step four: According to the maximum curvature radius of the head stagnation point of the sub-orbital aircraft, the minimum curvature radius of the head stagnation point, the curvature radius of the wing leading edge, the wing leading edge sweep angle, and the flight condition in step two, the head stagnation point heat flux density, the wing inner wing segment leading edge heat flux density, and the wing outer wing segment leading edge heat flux density of the sub-orbital aircraft are obtained.

[0049] The specific implementation method of step four is as follows:

[0050] ①. Based on equation (14) and equation (15), the stagnation point gas enthalpy and the gas wall enthalpy value are calculated

[0051]

[0052] In the equation, c p is the air specific heat ratio, T ∞ is the incoming flow temperature, R is the ideal gas constant, W is the atmospheric constant, and T w is the wall temperature.

[0053] ii. Calculate the stagnation heat flux value based on formula (16)

[0054]

[0055] In the formula, r min is the minimum head stagnation point curvature radius, r max is the maximum head stagnation point curvature radius, p sea is the sea level atmospheric density, p ∞ is the incoming flow density

[0056] iii. Calculate the inner and outer wing segment leading edge heat flux values based on formula (17), formula (18)

[0057]

[0058] In the formula, q sp is the equal radius sphere heat flux, x int is the wing inner wing segment leading edge sweepback angle, x ext is the wing outer wing segment leading edge sweepback angle, n a is the correction index, a is the angle of attack in the flight condition; q sp and n a The specific values are given by formula (19), formula (20)

[0059]

[0060] In the formula, r l is the wing leading edge curvature radius, h 300K is the air enthalpy value at a temperature of 300K, x is the wing sweepback angle

[0061] Step five: According to the sub-orbital vehicle overall volume, total surface area, wing area, and tail area in step two, obtain the reentry return mass of the sub-orbital vehicle.

[0062] The specific implementation method of step five is as follows:

[0063] i. Calculate the sub-orbital vehicle dry weight based on formula (21)

[0064] m dry = m ext + m p + m tps + m pl + m uc + m es + m of (21)

[0065] In the formula, m ext is the external structure mass, m p is the propulsion system mass, m tps is the thermal protection system mass, muc For landing gear mass, m es For the quality of onboard electronic systems, m of For other liquid masses, m ext Equation (22) represents

[0066]

[0067] In the formula, m fu For fuselage weight, m wing For wing mass, m vt For the tail fin mass, w fu For the fuselage unit mass factor, w wing For the wing's unit mass factor, w vt For the tail fin's unit mass factor, S fu For fuselage surface area, S wing For wing surface area, S vt Let g be the surface area of ​​the tail fin, and g0 be the gravitational acceleration at sea level.

[0068] Propulsion system mass m p Quality of thermal protection system m tps Fuel mass m fuel Other liquid mass m of Equations (23)-(26) represent

[0069] m p =n p ×w p ×S in (twenty three)

[0070] m tps =w tps ×S wing (twenty four)

[0071] m fuel =ρ fuel ×V tank ×k fuel (25)

[0072] m of =m fuel ×k of (26)

[0073] In the formula, n p For the number of engines, w p For engine unit mass factor, S in For the engine intake duct area, w tps For the unit mass factor of the thermal protection system, ρ fuel fuel density, V tank For fuel tank volume, k fuelfuel loading factor, k of fuel-other liquid ratio factor

[0074] ②. Solve the linear equations in equation (27) to obtain the suborbital vehicle takeoff mass m takeoff and dry mass m dry

[0075]

[0076] ③. Calculate the suborbital vehicle reentry return mass based on equation (28)

[0077] m entry = m dry +(f u +f a )·m fuel (28)

[0078] Step six: Establish the trajectory model of the suborbital vehicle according to the reference area of the suborbital vehicle in step two, the lift coefficient and the drag coefficient in step three, the suborbital vehicle reentry return mass, and the flight working condition.

[0079] The specific implementation method of step six is as follows:

[0080] ①. Determine the suborbital vehicle reentry return attack angle based on equation (29):

[0081]

[0082] In the formula, v(t) is the flight speed; a max is the initial reentry attack angle, a min is the final reentry attack angle, k α is the attack angle descent rate, v max is the maximum speed before the start of aerodynamic deceleration, t1 is the time required to reach the maximum speed, t2 is the time experienced when the speed change amount is Δv; Δv is represented by equation (30)

[0083]

[0084] ②. Solve the partial differential equations in equation (31) to obtain the speed, speed inclination angle, heading angle, latitude, longitude, and geocentric distance of the suborbital vehicle during the reentry return process

[0085]

[0086] In the formula, V is the flight speed, θ is the speed inclination angle, σ is the yaw angle, r is the geocentric distance, φ is the latitude, λ is the longitude, m is the reentry return mass of the vehicle, g r is the gravity acceleration geocentric direction component, and ω eωe is the angular velocity of the earth rotation, g ω is the component of the earth rotation direction of the gravity acceleration, v is the inclination angle of the suborbital vehicle, L(v, a) is the lift, and D(v, a) is the drag.

[0087] Step seven: based on the suborbital vehicle carrying capacity optimization model in step one, taking the maximum payload mass of the suborbital vehicle as the optimization objective, using the initial sample expansion mechanism to expand the suborbital vehicle sample, using all sample information in the sample database to construct the Kriging surrogate model of the objective function and the constraint function, using the Kriging surrogate model to replace the high-time-consuming analysis model of the suborbital vehicle, and solving the surrogate model pseudo-optimal solution by the constrained differential evolution algorithm; combining the comprehensive improvement probability and the key design space method, realizing the management and updating of the surrogate model under the condition of meeting the design requirements of each discipline of the suborbital vehicle, and efficiently guiding the convergence of the suborbital vehicle carrying capacity optimization process to the global optimal solution.

[0088] The specific implementation method of step seven is as follows:

[0089] ①. Determine the initial parameters of the approximate optimization strategy based on the initial sample expansion mechanism, including the initial sampling number n doe , the constraint threshold scaling factor η, the initial sample expansion number n l , the number of newly added sample points in the local search stage n a , and the maximum number of model calls

[0090] ②. Use the Latin hypercube experimental design method to construct n doe initial sample points in the design space, and call the true analysis model to calculate the model response values at the initial sample points; add all the sample points and response values to the sample database; set the optimization iteration number iter to 1;

[0091] ③. Based on the initial sample expansion mechanism, obtain n l sample points;

[0092] (1). Based on the sample database in step ②, establish the random forest classifier as shown in formula (32)

[0093]

[0094] In the formula, n RF is the number of sample points for constructing the random forest classifier, is the classification value of the i-th sample point, c i is the constraint violation degree of the i-th sample point, is the normalized value of the j-th constraint function, is the constraint threshold value of the j-th constraint function, and P thresh is the constraint violation threshold value; Pthresh is denoted as

[0095] P thresh is denoted as i,min + η · (c i,max - c i,min ) (33)

[0096] where c i,min is the minimum constraint violation, c i,max is the maximum constraint violation, and η is the constraint threshold scaling factor;

[0097] A large number of simple sample points are generated in the design space, and the simple sample points are classified by the constructed random forest classifier; if the classification value is 1, the simple sample point is recorded as a high-quality sample point;

[0098] (2). The number of high-quality sample points is counted; if there is no high-quality sample point, the constraint threshold is increased by 10%; otherwise, the high-quality sample point clustering center is obtained by the K-means method;

[0099] (3). A sample point is randomly generated between the high-quality sample point clustering center and the sample point with the minimum objective function response value in the sample database of step ②; the real analysis model is called to calculate the real model response value of the sample point, and the sample information is added to the sample database of step ②; if the number of new sample points reaches n l , the expansion mechanism is terminated; otherwise, step (1) is returned, and the expansion process continues;

[0100] ④. The Kriging surrogate model of the objective function and the constraint function is constructed using all the sample information in the sample database, and the constructed Kriging surrogate model is optimized by the constraint differential evolution method to obtain the pseudo-optimal solution The real analysis model is called to calculate the real model response value of the pseudo-optimal solution, and the sample information of the pseudo-optimal solution is added to the sample database;

[0101] ⑤. The sub-optimization problem shown in equation (34) is solved to obtain the sample point x PI

[0102]

[0103] where PI(·) is the objective function improvement probability, P(gi(·)≤0) is the i-th constraint satisfaction probability, d m (X,X P-min ) is the Manhattan distance between the sample point and the optimal comprehensive improvement probability sample point in the current sample database; PI(·), P(gi(·)≤0), and d m (X,X P-min) are represented by formula (35), formula (36), formula (37) respectively

[0104]

[0105] In the formula, Φ(·) is a distribution function of a standard normal distribution, f min is a minimum value of a target function in a sample database, is a Kriging proxy model prediction value of the target function, is a Kriging proxy model prediction value of the i-th constraint function, s(·) is a Kriging proxy model prediction variance of the target function, s g,i (·) is a Kriging proxy model prediction variance of the i-th constraint function; a sample point x PI with the maximum comprehensive improvement probability is calculated by calling a real analysis model, and x PI sample information is added to the sample database.

[0106] ⑥. Based on the current sample database, a key sampling space is constructed; the radius of the key sampling space is shown in formula (38)

[0107]

[0108] In the formula, x opt,k is the optimal feasible solution of the k-th iteration, x opt,k-1 is the optimal feasible solution of the k-1-th iteration, is a sample point with the maximum error obtained by the one-by-one checking method when the iteration number is 1; if the optimality of x opt,k is improved compared with x opt,k-1 , x opt,k is taken as the center of the key design space, otherwise x opt,k-1 is taken as the center of the key design space.

[0109] A sample point x SSS is randomly generated in the key sampling space, the real analysis model is called to calculate the real model response value of x SSS , and x SSS sample information is added to the sample database.

[0110] ⑦. It is judged whether the number of sample points in the current sample database reaches If it reaches, the optimization is stopped; otherwise, the iteration number of the optimization is iter+1, step seven ④ is turned to, the Kriging proxy models of the target function and the constraint function are updated, and the optimization process continues.

[0111] Step eight: the feasible optimal solution in the current sample database is output as the optimization result of the sub-orbital vehicle, that is, the sub-orbital vehicle carrying performance optimization considering flight state constraints is realized.

[0112] Further comprising step nine: the flight state constraint considering suborbital vehicle carrying performance optimization method described in steps one to eight is applied to suborbital vehicle carrying performance optimization, and the performance and design efficiency of the suborbital vehicle are improved; according to the suborbital vehicle optimization result obtained in step eight, the suborbital vehicle realizes the maximization of carrying performance under the constraints of take-off mass, maximum axial overload, maximum normal overload, re-entry segment terminal velocity, re-entry segment range, re-entry segment terminal lift-drag ratio, maximum flight dynamic pressure, re-entry segment terminal velocity inclination, maximum stagnation point heat flux density and maximum wing leading edge heat flux density. The suborbital vehicle carrying performance optimization includes single-stage horizontal take-off and landing suborbital vehicle carrying performance optimization, two-stage horizontal take-off and landing suborbital vehicle carrying performance optimization, vertical take-off and landing suborbital vehicle carrying performance optimization and other fields.

[0113] Beneficial effects:

[0114] 1. In view of the problem of time-consuming calculation of multi-disciplinary analysis model and insufficient data mining under limited calculation cost, the flight state constraint considering suborbital vehicle carrying performance optimization method disclosed in the present application comprehensively considers the parameter transmission relationship among the geometric main model, aerodynamic force and heat model, mass estimation model and trajectory model of the suborbital vehicle, and establishes a multi-disciplinary analysis model of the suborbital vehicle considering flight state constraints such as overload, dynamic pressure and heat flow. An approximate optimization strategy based on initial sample expansion mechanism is adopted to realize efficient mining of high-quality design schemes of the suborbital vehicle. The flight state constraint considering suborbital vehicle carrying performance optimization method disclosed in the present application combines an approximate optimization strategy based on initial sample expansion mechanism, expands the design variables of the suborbital vehicle, constructs Kriging surrogate models of the objective function and constraint function by using all sample information in the sample database, uses the Kriging surrogate model to replace the high-time-consuming analysis model of the suborbital vehicle, and solves the pseudo-optimal solution of the surrogate model by a constraint differential evolution algorithm; combined with the comprehensive improvement probability and key design space method, the surrogate model management and updating are realized under the condition of meeting the design requirements of each discipline of the suborbital vehicle, and the suborbital vehicle carrying performance optimization process is efficiently guided to converge to the global optimal solution. The effective load of the optimization scheme obtained by the method can be improved by more than 171 kg compared with the initial scheme.

[0115] 2. The suborbital vehicle carrying performance optimization method considering flight state constraints, which adopts an approximate optimization strategy based on an initial sample expansion mechanism to expand the suborbital vehicle design variables, uses all sample information in a sample database to construct a Kriging surrogate model of a target function and a constraint function, adopts the Kriging surrogate model to replace a high-time-consuming analysis model of the suborbital vehicle, and solves a pseudo-optimal solution of the surrogate model by a constraint differential evolution algorithm; in combination with a comprehensive improvement probability and a key design space method, the surrogate model management and update are realized under the condition that the design requirements of each discipline of the suborbital vehicle are met, the suborbital vehicle carrying performance optimization process is efficiently guided to converge to a global optimal solution, and the efficiency of the suborbital vehicle carrying performance optimization is further improved. The global optimal solution corresponds to a suborbital vehicle optimization result, and the suborbital vehicle realizes the maximization of carrying performance under the condition of flight state constraints according to the suborbital vehicle optimization result. BRIEF DESCRIPTION OF DRAWINGS

[0116] Figure 1 It is a suborbital vehicle carrying performance optimization method considering flight state constraints according to the present application;

[0117] Figure 2 It is a suborbital vehicle geometric configuration;

[0118] Figure 3 It is a flowchart of an approximate optimization strategy based on an initial sample expansion mechanism;

[0119] Figure 4 It is a suborbital vehicle optimization convergence curve diagram;

[0120] Figure 5 It is a comparison of suborbital vehicle configurations before and after optimization;

[0121] Figure 6 It is a comparison of suborbital vehicle flight states, wherein Figure 6 (a) is a suborbital vehicle altitude-time curve diagram, Figure 6 (b) is a suborbital vehicle flight Mach number curve diagram, Figure 6 (c) is a suborbital vehicle axial overload curve diagram, Figure 6 (d) is a suborbital vehicle normal overload curve diagram, Figure 6 (e) is a suborbital vehicle flight dynamic pressure curve diagram, Figure 6 (f) is a suborbital vehicle stagnation point heat flux density curve diagram, Figure 6 (g) is a suborbital vehicle wing leading edge heat flux density curve diagram. DETAILED DESCRIPTION

[0122] In order to further illustrate the purposes and advantages of the present application, the present application will be further described below in combination with the drawings and examples.

[0123] Embodiment 1: Commercial suborbital launch vehicle multidisciplinary design optimization example

[0124] The suborbital launch vehicle carrying performance optimization method considering flight state constraints disclosed in this embodiment is suitable for commercial suborbital launch vehicle multidisciplinary design optimization problems, can guarantee obtaining a high-performance design scheme of a suborbital launch vehicle in the overall design stage, and shortens the design cycle of the suborbital launch vehicle.

[0125] The suborbital launch vehicle carrying performance optimization method considering flight state constraints disclosed in this embodiment specifically implements the following steps:

[0126] Step one: determining a suborbital launch vehicle carrying performance optimization model and initial parameters of an approximate optimization strategy based on an initial sample expansion mechanism.

[0127] The specific implementation method of step one is as follows:

[0128] Taking the suborbital launch vehicle payload mass m pl as the objective function, and considering the takeoff mass m takeoff , axial overload N a , normal overload N n , terminal velocity of reentry segment V end , reentry segment range R, terminal lift-drag ratio C l / C d of reentry segment, flight dynamic pressure Q, terminal velocity inclination angle Θ touch of reentry segment, stagnation point heat flux density q stag , and wing leading edge heat flux density q wing constraints, a suborbital launch vehicle carrying performance optimization model of the suborbital launch vehicle is established, as shown in formula (39):

[0129] find X=[c v,1 ,c v,2 ,c h,1 ,c h,2 ,L int ,L ext ,b root ,b tip ,χ int ,χ ext ,L h ,W r ,f u ,f a ,k α ]

[0130] min f(X)=-m pl

[0131]

[0132] In the formula, X represents a design variable, specifically including aircraft head profile vertical direction control point 1 coordinate c v,1 , aircraft head profile vertical direction control point 2 coordinate c v,2 , aircraft head profile horizontal direction control point 1 coordinate c h,1 , aircraft head profile vertical direction control point 2 coordinate c h,2 , inner wing segment half span length L int , outer wing segment half span length L ext , wing root chord length b root , wing tip chord length b tip , inner wing segment sweepback angle χ int , outer wing segment sweepback angle χ ext , head length L h , afterbody width W r , unpredictable fuel mass ratio f u , additional fuel mass ratio f a , and attack angle reduction rate k α as design variables. X LB and X UB respectively represent lower and upper limits of the design variables.

[0133] The initial parameters of the approximate optimization strategy based on the initial sample expansion mechanism include initial sampling number n doe , constraint threshold scaling factor η, initial sample expansion number n l , number of newly added sample points in the local search stage n a , and maximum model call number

[0134] Step two: based on the design variables of the suborbital aircraft, a parameterized geometric main model of the suborbital aircraft is established, and the geometric characteristics (reference area, head stagnation point maximum curvature radius, head stagnation point minimum curvature radius, wing leading edge curvature radius, wing leading edge sweepback angle, total volume, total surface area, wing area, and tail area) of the suborbital aircraft are obtained.

[0135] The specific implementation method of step two is as follows:

[0136] Based on the NURBS method shown in formula (40), the head, fuselage, wing, afterbody, and tail of the suborbital aircraft are parameterized modeled to obtain the parameterized geometric main model:

[0137]

[0138] In the formula, d i is a curve control point, R i,p (u) is a rational basis function defined on u∈[0,1], ω i is a weight coefficient, and N i,p(u) is the i-th p-th order B-spline basis function, n represents the number of nodes. N i,p (u) can be represented by formula (41), formula (42); when the node u and the order p are determined, the B-spline basis function can be uniquely determined

[0139]

[0140] In the formula, N i,p (u) is the i-th p-th order B-spline basis function, u i is the i-th node.

[0141] The parameterized shapes of the head, the body, the wing, the rear body and the tail are generated by formula (40)-(42); and the geometric characteristics (reference area, maximum curvature radius of head stagnation point, minimum curvature radius of head stagnation point, wing leading edge curvature radius, wing leading edge sweep angle, total volume, total surface area, wing area and tail area) are measured as inputs for obtaining the lift coefficient, the drag coefficient, the head stagnation point heat flux density, the wing inner wing segment leading edge heat flux density, the wing outer wing segment leading edge heat flux density, the reentry return mass and the trajectory model.

[0142] Step three: according to the sub-orbital vehicle geometric master model and the reference area in step two, the lift coefficient and the drag coefficient of the sub-orbital vehicle are obtained.

[0143] The specific implementation method of step three is as follows:

[0144] ①. Based on the parameterized geometric master model and the calculation domain in step two, unstructured grids are divided, and the grid size is 0.5m; and local refinement processing is added at the head stagnation point and the wing, and the grid size is 0.15m.

[0145] Boundary layer grids are added at the head and the wing, and the first layer grid height is shown in formula (43)

[0146]

[0147] In the formula, the dimensionless wall distance y + ≈1, L r is the reference length of the vehicle, Re is the Reynolds number, which is represented by formula (44)

[0148]

[0149] In the formula, p ∞ is the atmospheric density, v ∞ is the incoming flow velocity in the flight working condition, and μ is the atmospheric viscosity.

[0150] After the local refinement processing and the addition of the boundary layer grids, the grid A is obtained.

[0151] ②. Based on the reference area of the sub-orbital vehicle measured in step two and the grid A in ①, the lift coefficient and the drag coefficient of the sub-orbital vehicle are solved by the CFD method. In the CFD solution, the density-based solver and the energy equation are used, and the far-field pressure boundary condition is set for the incoming flow; the convection upwind splitting format is used for spatial discretization, and the Green-Gauss node method is used to calculate the gradient; in addition, 3-level FMG is used for initialization before starting the solution to improve the solution convergence speed.

[0152] Step four: According to the maximum curvature radius of the head point of the sub-orbital vehicle, the minimum curvature radius of the head point of the sub-orbital vehicle, the curvature radius of the wing leading edge, the sweep angle of the wing leading edge, and the flight working condition in step two, the heat flux density of the head point of the sub-orbital vehicle, the heat flux density of the leading edge of the inner wing section of the wing, and the heat flux density of the leading edge of the outer wing section of the wing are obtained.

[0153] The specific implementation method of step four is as follows:

[0154] ①. Based on formula (45) and formula (46), the stagnation point gas enthalpy and the gas wall enthalpy value are calculated

[0155]

[0156] In the formula, c p is the specific heat ratio of air, T ∞ is the incoming flow temperature, R is the ideal gas constant, W is the atmospheric constant, T w is the wall temperature;

[0157] ②. Based on formula (47), the stagnation point heat flux density value is calculated

[0158]

[0159] In the formula, r min is the minimum head point curvature radius, r max is the maximum head point curvature radius, p sea is the sea level atmospheric density, p ∞ is the incoming flow density;

[0160] ③. Based on formula (48) and formula (49), the inner and outer wing section leading edge heat flux density values are calculated

[0161]

[0162]

[0163] In the formula, q sp is the equal radius ball heat flux, x int is the wing inner wing section leading edge sweep angle, x ext is the wing outer wing section leading edge sweep angle, n awherein a is the angle of attack in flight condition, q is the dynamic pressure sp wherein n a The specific value of n

[0164]

[0165] wherein r l is the wing leading edge radius, h 300K is the air enthalpy at 300 K, and χ is the wing sweep angle.

[0166] Step five: Obtain the reentry return mass of the suborbital vehicle according to the total volume, total surface area, wing area, and tail area of the suborbital vehicle in step two.

[0167] The specific implementation method of step five is as follows:

[0168] ①. Calculate the dry mass of the suborbital vehicle based on formula (52)

[0169] m dry = m ext + m p + m tps + m pl + m uc + m es + m of (52)

[0170] wherein m ext is the external structure mass, m p is the propulsion system mass, m tps is the thermal protection system mass, m uc is the landing gear mass, m es is the on-board electronic system mass, and m of is the other liquid mass. Among them, m ext is represented by formula (53)

[0171]

[0172] wherein m fu is the fuselage mass, m wing is the wing mass, m vt is the tail mass, w fu is the fuselage unit mass factor, w wing is the wing unit mass factor, w vt is the tail unit mass factor, S fu is the fuselage surface area, S wing is the wing surface area, S vt is the tail surface area, and g0 is the sea level gravitational acceleration.

[0173] The propulsion system mass mp , the heat protection system mass m tps , the fuel mass m fuel and other liquid mass m of are expressed by formula (54) - formula (57)

[0174] m p = n p × w p × S in (54)

[0175] m tps = w tps × S wing (55)

[0176] m fuel = ρ fuel × V tank × k fuel (56)

[0177] m of = m fuel × k of (57)

[0178] In the formula, n p is the engine number, w p is the engine unit mass factor, S in is the engine air inlet area, w tps is the heat protection system unit mass factor, ρ fuel is the fuel density, V tank is the fuel tank volume, k fuel is the fuel loading coefficient, k of is the fuel-other liquid ratio coefficient;

[0179] ②. Solve the linear equations in formula (58) to obtain the suborbital vehicle takeoff mass m takeoff and the dry mass m dry

[0180]

[0181] ③. Calculate the suborbital vehicle reentry return mass based on formula (59)

[0182] m entry = m dry + (f u + f a )·m fuel (59)

[0183] Step six: according to the suborbital vehicle reference area in step two, the lift coefficient and the drag coefficient in step three, the suborbital vehicle reentry return mass and the flight condition, the trajectory model of the suborbital vehicle is established.

[0184] The specific implementation method of step six is as follows:

[0185] ①. Based on formula (60), the suborbital vehicle reentry return attack angle is determined:

[0186]

[0187] In the formula, v(t) is the flight speed; α max is the initial reentry attack angle, α min is the final reentry attack angle, k α is the attack angle drop rate, v max is the maximum speed before the start of aerodynamic deceleration, t1 is the time required to reach the maximum speed, t2 is the time experienced when the speed change is Δv; Δv is represented by formula (61)

[0188]

[0189] ②. The partial differential equation system in formula (62) is solved to obtain the speed, speed inclination angle, heading angle, latitude, longitude and geocentric distance of the suborbital vehicle during reentry return

[0190]

[0191] In the formula, V is the flight speed, θ is the speed inclination angle, σ is the yaw angle, r is the geocentric distance, φ is the latitude, λ is the longitude, m is the reentry return mass of the vehicle, g r is the geocentric component of gravitational acceleration, ω e is the earth rotation angular velocity, g ω is the geodetic component of gravitational acceleration, v is the suborbital vehicle roll angle, L(v, α) is the lift, and D(v, α) is the drag.

[0192] Step seven: based on the suborbital vehicle carrying performance optimization model in step one, the maximum payload mass of the suborbital vehicle is taken as the optimization objective, and the approximate optimization strategy based on the initial sample expansion mechanism is used to optimize the suborbital vehicle design variables. Under the condition of meeting the requirements of various disciplines of the suborbital vehicle design, the efficient exploration of the maximum payload scheme is realized.

[0193] The specific implementation method of step seven is as follows:

[0194] ①. Determine the initial parameters of the approximate optimization strategy based on the initial sample expansion mechanism, including the initial sampling number n doe , the constraint threshold scaling factor η, and the initial sample expansion number nl , the number of newly added sample points in the local search stage a , and the maximum number of model calls

[0195] ②. Construct n doe initial sample points in the design space using the Latin hypercube design method, and call the real analysis model to calculate the model response values at the initial sample points; add all sample points and response values to the sample database; set the optimization iteration number iter to 1;

[0196] ③. Obtain n l sample points based on the initial sample expansion mechanism;

[0197] (1). Based on the sample database in step ②, establish a random forest classifier as shown in formula (63)

[0198]

[0199] In the formula, n RF is the number of sample points for constructing the random forest classifier, is the classification value of the i-th sample point, c i is the constraint violation degree of the i-th sample point, is the normalized value of the j-th constraint function, is the constraint threshold of the j-th constraint function, P thresh is the constraint violation threshold; P thresh represents

[0200] P thresh = c i,min + η · (c i,max - c i,min ) (64)

[0201] In the formula, c i,min is the minimum constraint violation degree, c i,max is the maximum constraint violation degree, and η is the constraint threshold scaling factor.

[0202] Generate 10,000 simple sample points in the design space, and classify the simple sample points by the constructed random forest classifier; if the classification value is 1, record the simple sample point as a high-quality sample point;

[0203] (2). Count the number of high-quality sample points; if there are no high-quality sample points, increase the constraint threshold by 10%; otherwise, obtain the high-quality sample point clustering center by the K-means method;

[0204] (3). Randomly generate a sample point between the clustering center of high-quality sample points and the sample point with the minimum objective function response value in the sample database of step ②; call the real analysis model to calculate the real model response value of the sample point, and add the sample information of the sample point to the sample database of step ②; check whether the number of newly added sample points reaches n l , if yes, the expansion mechanism terminates; otherwise, return to step (1) and continue the expansion process;

[0205] ④. Use all sample information in the sample database to construct the Kriging surrogate model of the objective function and the constraint function, and optimize the constructed Kriging surrogate model by the constrained differential evolution method to obtain a pseudo-optimal solution Call the real analysis model to calculate the real model response value of the pseudo-optimal solution, and add the sample information of the pseudo-optimal solution to the sample database;

[0206] ⑤. Solve the sub-optimization problem shown in formula (65) to obtain a sample point x PI

[0207]

[0208] In the formula, PI(·) is the objective function improvement probability, P(gi(·)≤0) is the i-th constraint satisfaction probability, d m (X,X P-min ) is the Manhattan distance between the sample point and the optimal comprehensive improvement probability sample point in the current sample database; PI(·), P(gi(·)≤0) and d m (X,X P-min ) are represented by formula (66), formula (67) and formula (68) respectively

[0209]

[0210] In the formula, Φ(·) is the distribution function of the standard normal distribution, f min is the minimum value of the objective function in the sample database, is the Kriging surrogate model prediction value of the objective function, is the Kriging surrogate model prediction value of the i-th constraint function, s(·) is the Kriging surrogate model prediction variance of the objective function, s g,i (·) is the Kriging surrogate model prediction variance of the i-th constraint function; call the real analysis model to calculate the real model response value of the sample point with the maximum comprehensive improvement probability, and add the sample information of the sample point to the sample database; PI ; PI

[0211] ​⑥. Based on the current sample database, a key sampling space is constructed; the radius of the key sampling space is shown as formula (69)

[0212]

[0213] In the formula, x opt,k is the optimal feasible solution of the kth iteration, x opt,k-1 is the optimal feasible solution of the (k-1)th iteration, is the sample point with the largest error obtained by the one-by-one checking method when the iteration number is 1; if the optimality of x opt,k is improved compared with x opt,k-1 , x opt,k is taken as the center of the key design space, otherwise x opt,k-1 is taken as the center of the key design space

[0214] A sample point x SSS is randomly generated in the key sampling space, the real analysis model is called to calculate the real model response value of x SSS , and the sample information of x SSS is added to the sample database

[0215] ⑦. It is judged whether the number of sample points in the current sample database reaches If yes, the optimization is stopped; otherwise, the iteration number iter of the optimization is iter+1, step seven (4) is turned to, the Kriging surrogate model of the objective function and the constraint function is updated, and the optimization process continues.

[0216] Step eight: the feasible optimal solution in the current sample database is output as the optimization result of the suborbital vehicle design scheme.

[0217] In order to better reflect the effectiveness and engineering practicability of the present application, the following takes the multi-disciplinary design optimization problem of a commercial suborbital vehicle as an example to further illustrate the present application in combination with the drawings and tables.

[0218] Step nine: in the embodiment, the value range of the design variable is c v,1 ∈[750mm,900mm], c v,2 ∈[2000mm,2200mm], c h,1 ∈[1150mm,1350mm], c h,2 ∈[2100mm,2300mm], L int ∈[2200mm,5000mm], L ext ∈[3800mm,6500mm], b root ∈[10000mm,12300mm], b tip ∈[1200mm,1800mm], χ int∈[50°,60°], χ ext ∈[15°,25°], L h ∈[4230mm,5170mm], W r ∈[1980mm,2420mm], f u ∈[0.018,0.022], f a ∈[0.018,0.022], k α ∈[0.010°s / m,0.015°s / m]. The approximate optimization strategy based on the initial sample expansion mechanism is set as follows: the initial sampling number n doe is 30, the constraint threshold scaling factor η is 0.2, the initial sample expansion number n l is 30, the number of newly added sample points in the local search stage n a is 1, and the maximum model calling number is 150. The optimization convergence curve obtained by the approximate optimization strategy based on the initial sample expansion mechanism and the efficient global optimization method is shown in FIG. 3, the target function comparison before and after optimization is shown in Table 1, the design variables are shown in Table 2, and the constraint conditions are shown in Table 3 Figure 4

[0219] Table 1. Target function comparison before and after optimization

[0220]

[0221] Table 2. Design variable comparison before and after optimization

[0222]

[0223] Table 3. Constraint condition comparison before and after optimization

[0224]

[0225]

[0226] As shown in Table 1, the payload mass of the optimization scheme is increased by more than 171 kg compared with the initial scheme. As shown in Table 2 and Figure 5 , the configuration of the optimization scheme changes greatly compared with the configuration of the initial scheme. The span length and wing root chord length are increased by 789.07 mm and 1054.89 mm, respectively, thereby providing more lift. As shown in Figure 6 (b), the Mach number of the optimization scheme in the reentry return stage is lower than that of the initial scheme, thereby having higher lift and lift-drag ratio. As shown in Figure 6 (e), the maximum dynamic pressure of the optimization scheme is reduced by 23.21% compared with the initial scheme, thereby reducing the structural design burden of the suborbital spacecraft. In Figure 6 (f) and Figure 6 ​​In the (g), the optimization scheme reduces the maximum heat flow at the stagnation point and the leading edge of the wing by 64.85% and 14.95% respectively compared with the initial scheme, thereby reducing the ablation risk faced by the suborbital vehicle in the reentry return process. In summary, the above optimization results show that the suborbital vehicle carrying performance optimization method considering flight state constraints disclosed in the present application can obtain a suborbital vehicle design scheme with higher payload and meeting the actual engineering requirements compared with the initial scheme. The above optimization results verify the rationality, effectiveness and engineering practicability of the present application.

[0227] The specific description described above further details the purpose, technical scheme and beneficial effects of the application. It should be understood that the above description is only a specific embodiment of the application and is not used to limit the protection scope of the application. Any modification, equivalent replacement, improvement, etc. made within the spirit and principles of the application shall be included in the protection scope of the application.

Claims

1. A method for suborbital vehicle performance optimization considering flight state constraints, characterized in that: Comprising the following steps, Step one: determine the sub-orbital vehicle carrying performance optimization model, and determine the initial parameters of the approximate optimization strategy based on the initial sample expansion mechanism; Step two: based on the sub-orbital vehicle design variables, establish the parameterized geometric main model of the sub-orbital vehicle, and obtain the geometric characteristics of the sub-orbital vehicle; The geometric characteristics of the sub-orbital vehicle include reference area, maximum curvature radius of head point, minimum curvature radius of head point, wing leading edge curvature radius, wing leading edge sweep angle, total volume, total surface area, wing area and tail area; Step three: according to the geometric main model of the sub-orbital vehicle in step two and the reference area, obtain the lift coefficient and drag coefficient of the sub-orbital vehicle; Step four: according to the maximum curvature radius of the head point, the minimum curvature radius of the head point, the wing leading edge curvature radius, the wing leading edge sweep angle of the sub-orbital vehicle in step two, and the flight working condition, obtain the head point heat flux density, the wing inner wing segment leading edge heat flux density and the wing outer wing segment leading edge heat flux density of the sub-orbital vehicle; Step five: according to the total volume, total surface area, wing area and tail area of the sub-orbital vehicle in step two, obtain the reentry return mass of the sub-orbital vehicle; Step six: according to the reference area of the sub-orbital vehicle in step two, the lift coefficient and the drag coefficient in step three, the reentry return mass of the sub-orbital vehicle and the flight working condition, establish the trajectory model of the sub-orbital vehicle; Step seven: based on the sub-orbital vehicle carrying performance optimization model in step one, taking the maximum payload mass of the sub-orbital vehicle as the optimization objective, using the initial sample expansion mechanism to expand the sub-orbital vehicle sample, using all sample information in the sample database to construct the Kriging surrogate model of the objective function and the constraint function, using the Kriging surrogate model to replace the high time-consuming analysis model of the sub-orbital vehicle, and solving the pseudo optimal solution of the surrogate model by the constrained differential evolution algorithm; Combined with the comprehensive improvement probability and the key design space method, the proxy model management and update are realized under the condition of meeting the design requirements of each discipline of the sub-orbital vehicle, and the sub-orbital vehicle carrying performance optimization process is efficiently guided to converge to the global optimal solution; Step eight: output the feasible optimal solution in the current sample database as the optimization result of the sub-orbital vehicle, that is, realize the sub-orbital vehicle carrying performance optimization considering flight state constraints.

2. The method of claim 1, wherein: The specific implementation method of step one is as follows: Determine the sub-orbital vehicle sub-orbital vehicle carrying performance optimization model and the initial parameters of the approximate optimization strategy based on the initial sample expansion mechanism; Maximum target function, taking into account the take-off mass m pl Maximum target function, taking into account the take-off mass m takeoff Maximum target function, taking into account the take-off mass m a Maximum target function, taking into account the take-off mass m n Maximum target function, taking into account the take-off mass m end Maximum target function, taking into account the take-off mass m l Maximum target function, taking into account the take-off mass m d Maximum target function, taking into account the take-off mass m touch Maximum target function, taking into account the take-off mass m stag Maximum target function, taking into account the take-off mass m wing The suborbital vehicle carrying performance optimization model of the suborbital vehicle is established according to the constraints, as shown in formula (1): find X=[c v,1 ,c v,2 ,c h,1 ,c h,2 ,L int ,L ext ,b root ,b tip ,χ int ,χ ext ,L h ,W r ,f u ,f a ,k α ] min f(X) = -m pl In the formula, X represents a design variable, and specifically includes a vertical direction control point 1 coordinate c of an aircraft head profile v,1 a vertical direction control point 2 coordinate c of the aircraft head profile v,2 a horizontal direction control point 1 coordinate c of the aircraft head profile h,1 a vertical direction control point 2 coordinate c of the aircraft head profile h,2 an inner wing segment half-span length L int an outer wing segment half-span length L ext a wing root chord length b root a wing tip chord length b tip an inner wing segment sweepback angle χ int an outer wing segment sweepback angle χ ext a head length L h a body width W r an unexpected fuel mass ratio f u an additional fuel mass ratio f a and an angle of attack descent rate k α X represents a design variable; X LB and X UB respectively represent a lower limit and an upper limit of the design variable. The initial parameter of the approximate optimization strategy based on the initial sample expansion mechanism includes an initial sampling number n doe , a constraint threshold scaling factor η, an initial sample expansion number n l , a number of newly added sample points in the local search stage n a , and a maximum model call number 3. The method of claim 1, wherein: The specific implementation method of step two is as follows: Based on the NURBS method shown in formula (2), the head, fuselage, wing, afterbody and tail of the sub-orbital vehicle are parameterized modeled to obtain the parameterized geometric main model: where d i is a curve control point, R i,p (u) is a rational basis function defined on u∈[0,1], ω i is a weight coefficient, N i,p (u) is the i-th p-th order B-spline basis function, n represents the number of nodes; N i,p (u) can be represented by equation (3), equation (4); when the node u and the order p are determined, the B-spline basis function can be uniquely determined where N i,p (u) is the i-th p-th order B-spline basis function, u i is the i-th node; Parametric shapes of the head, fuselage, wings, afterbody and tail are generated from the equations (2)-(4); geometric features are measured, which are used as inputs to obtain the lift coefficient, drag coefficient, head stagnation point heat flux density, wing inboard section leading edge heat flux density, wing outboard section leading edge heat flux density, reentry return mass and trajectory model; the geometric features include reference area, head stagnation point maximum radius of curvature, head stagnation point minimum radius of curvature, wing leading edge radius of curvature, wing leading edge sweep angle, total volume, total surface area, wing area and tail area.

4. The method of claim 1, wherein: The specific implementation method of step three is as follows: ①. Based on the parametric geometric main model and the calculation domain in step two, unstructured grids are divided; and local refinement processing is added at the head stagnation point and the wing; Boundary layer grids are added at the head and the wing, and the height of the first layer of grid is shown in equation (5) where y is a dimensionless wall distance + ≈1, L r L is a reference length for the aircraft and Re is the Reynolds number, given by equation (6) In the formula, ρ ∞ atmospheric density, v ∞ The velocity of the incoming flow during flight is μ, and the atmospheric viscosity is μ. After local refinement processing and adding boundary layer grids, grid A is obtained; ②. Based on the reference area of the suborbital vehicle measured in step two and the grid A in ①, the lift coefficient and the drag coefficient of the suborbital vehicle are solved by the CFD method; in the CFD solution, the density-based solver and the energy equation are used, and the incoming flow is set as the pressure far-field boundary condition; the convective upwind splitting format is used for spatial discretization, and the Green-Gauss node method is used to calculate the gradient; 3-level FMG is used for initialization before starting the solution to improve the solution convergence speed.

5. The method of claim 1, wherein: The specific implementation method of step four is as follows: ①. Based on equations (7) and (8), the stagnation point gas sensible enthalpy and the gas wall enthalpy value are calculated where c p is the air specific heat ratio, T ∞ is the incoming flow temperature, R is the ideal gas constant, W is the atmospheric constant, T w is the wall temperature; ②. Based on equation (9), the stagnation point heat flux density value is calculated where r min is the minimum head point radius of curvature, r max is the maximum head point radius of curvature, p sea is the sea level atmospheric density, p ∞ is the incoming flow density; ③. Based on equations (10) and (11), the inboard and outboard section leading edge heat flux density values are calculated where q sp is the constant heat flux, χ int is the sweep angle of the inner wing section, χ ext is the sweep angle of the outer wing section, n a is the correction exponent, α is the angle of attack in flight condition; q sp and n a are given by equations (12) and (13) respectively where r l is the wing leading edge radius of curvature, h 300K is the air enthalpy at 300 K, and χ is the wing sweep angle.

6. The method of claim 1, wherein: The specific implementation method of step five is as follows: ①. Based on equation (14), the dry weight of the suborbital vehicle is calculated m dry = m ext + m p + m tps + m pl + m uc + m es + m of (14) where m ext is the external structure mass, m p is the propulsion system mass, m tps is the thermal protection system mass, m uc is the landing gear mass, m es is the on-board electronic system mass, m of is the other liquids mass; where m ext is represented by equation (15) where m fu is the fuselage mass, m wing is the wing mass, m vt is the tail mass, w fu is the fuselage mass per unit factor, w wing is the wing mass per unit factor, w vt is the tail mass per unit factor, S fu is the fuselage surface area, S wing is the wing surface area, S vt is the tail surface area, g0is the sea level gravitational acceleration; Propulsion system mass m p Thermal protection system mass m tps Fuel mass m fuel And other liquid mass m of Is expressed by formula (16) - formula (19) m p = n p x w p x S in (16) m tps = w tps x S wing (17) m fuel = p fuel x V tank x k fuel (18) m of = m fuel x k of (19) where n p is the number of engines, w p is the engine unit mass factor, S in is the engine intake area, w tps is the heat shield unit mass factor, p fuel is the fuel density, V tank is the fuel tank volume, k fuel is the fuel loading factor, k of is the fuel-other liquid ratio factor; ii. Solving the linear equations in equation (20) to obtain the suborbital launch vehicle takeoff mass m takeoff and the dry mass m dry ③. Based on equation (21), the reentry return mass of the suborbital vehicle is calculated m entry = m dry + (f u + f a ) · m fuel (21 ) 。 7. The method of claim 1, wherein: The specific implementation method of step six is as follows: ①. Based on equation (22), the reentry return attack angle of the suborbital vehicle is determined: where v(t) is the flight velocity; a max is the initial reentry angle of attack, a min is the final reentry angle of attack, k α is the angle of attack rate, v max is the maximum velocity before aerodynamic deceleration begins, t1 is the time required to reach the maximum velocity, and t2 is the time experienced for a velocity change of Δv. Δv is represented by equation (23) ②. The partial differential equation system in equation (24) is solved to obtain the speed, speed inclination angle, heading angle, latitude, longitude and geocentric distance of the suborbital vehicle during the reentry return process In the formula, V is the flight speed, θ is the velocity tilt angle, σ is the yaw angle, r is the distance from the Earth's center, φ is the latitude, λ is the longitude, m is the reentry mass of the spacecraft, and g is the reentry return mass. r The component of gravitational acceleration in the geocentric direction, ω e For the Earth's rotational angular velocity, g ω Let ν be the gravitational acceleration component in the direction of Earth's rotation, ν be the suborbital tilt angle, L(v,α) be the lift, and D(v,α) be the drag.

8. The method of claim 1, wherein: The specific implementation method of step seven is as follows: ①. Determine the initial parameters of the approximate optimization strategy based on the initial sample expansion mechanism, including the initial sample number n doe , the constraint threshold scaling factor η, the initial sample expansion number n l , the number of newly added sample points in the local search stage n a , and the maximum number of model calls ②. Construct n within the design space using the Latin hypersquare experimental design method. doe 1. Establish initial sample points and call the real analysis model to calculate the model response value at the initial sample points; add all sample points and response values ​​to the sample database; The optimization iteration number iter is set to 1; ③. Based on the initial sample expansion mechanism, obtain n l sample points; (1). Based on the sample database in step ②, a random forest classifier as shown in equation (25) is established where n RF the number of sample points for constructing the random forest classifier, the classification value of the i-th sample point, c i the constraint violation degree of the i-th sample point, the normalized value of the j-th constraint function, the constraint threshold of the j-th constraint function, P thresh the constraint violation threshold; P thresh is represented as P thresh = c i,min + η · (c i,max - c i,min ) (26) where c i,min is the minimum constraint violation, c i,max is the maximum constraint violation, and η is a constraint threshold scaling factor. A large number of simple sample points are generated in the design space, and the simple sample points are classified by the established random forest classifier; if the classification value is 1, the simple sample point is recorded as a high-quality sample point; (2). The number of high-quality sample points is counted; if there is no high-quality sample point, the constraint threshold is increased by 10%; otherwise, the high-quality sample point clustering center is obtained by the K-means method; (3). In the high-quality sample point clustering center and the sample point with the minimum objective function response value in the sample database of step ②, a sample point is randomly generated; the real analysis model is called to calculate the real model response value of the sample point, and the sample information is added to the sample database of step ②; whether the number of newly added sample points reaches n is checked l , if yes, the expansion mechanism terminates; otherwise, step (1) is returned, and the expansion process continues; ④.Using all sample information in the sample database to construct the Kriging surrogate model of objective function and constraint function, and using the constrained differential evolution method to optimize the constructed Kriging surrogate model to obtain the pseudo-optimal solution The real analysis model is called to calculate the real model response value of the pseudo-optimal solution, and the sample information of the pseudo-optimal solution is added to the sample database. • Solve the sub-optimization problem shown in equation (27) to obtain the sample point x that maximizes the overall improvement probability within the design space PI find X=[c v,1 ,c v,2 ,c h,1 ,c h,2 ,L int ,L ext ,b root ,b tip ,χ int ,χ ext ,L h ,W r ,f u ,f a ,k α ] s.t.X lb ≤X≤X ub wherein PI(·) is an objective function improvement probability, P(g i (·)≤0) is an ith constraint satisfaction probability, d m (X,X P-min ) is a Manhattan distance between a sample point and an optimal comprehensive improvement probability sample point in a current sample database; PI(·), P(g i (·)≤0), and d m (X,X P-min ) are represented by equations (28), (29), and (30), respectively where Φ(·) is the distribution function of the standard normal distribution, f min is the minimum value of the objective function in the sample database, is the Kriging surrogate model prediction value of the objective function, is the Kriging surrogate model prediction value of the ith constraint function, s(·) is the Kriging surrogate model prediction variance of the objective function, s g,i (·) is the Kriging surrogate model prediction variance of the ith constraint function; the sample point x PI with the maximum comprehensive improvement probability is calculated by calling the real analysis model, and the real model response value of x PI is obtained; and the sample information is added to the sample database. ⑥. Based on the current sample database, a key sampling space is constructed; the radius of the key sampling space is shown in equation (31) where x opt,k is the optimal feasible solution of the kth iteration, x opt,k-1 is the optimal feasible solution of the (k-1)th iteration, is the sample point with the largest error obtained by the one-by-one checking method when the iteration number is 1; if the optimality of x opt,k is improved compared to x opt,k-1 , x opt,k is taken as the center of the key design space, otherwise x opt,k-1 is taken as the center of the key design space; Generate a sample point x randomly in the importance sampling space SSS , call the real analysis model to calculate the response value of x SSS , and add the sample information of x SSS to the sample database; ⑦. Determine whether the number of sample points in the current sample database reaches If yes, stop the optimization; otherwise, let the optimization iteration number iter be iter+1, go to step seven ④, update the Kriging surrogate model of the objective function and the constraint function, and continue the optimization process.

9. The method of claim 8, wherein: The objective function Kriging surrogate model in step seven is: payload mass m pl of the Kriging surrogate model; the constraint function Kriging surrogate model in step seven is: takeoff mass m takeoff of the Kriging surrogate model, axial overload N a of the Kriging surrogate model, normal overload N n of the Kriging surrogate model, terminal velocity V end of the Kriging surrogate model, terminal range R of the Kriging surrogate model, terminal lift-drag ratio C l / C d of the Kriging surrogate model, flight dynamic pressure Q of the Kriging surrogate model, terminal velocity inclination Θ touch of the Kriging surrogate model, nose stagnation point heat flux density q stag of the Kriging surrogate model, and wing leading edge heat flux density q wing of the Kriging surrogate model.

10. The method of claim 1, 2, 3, 4, 5, 6, 7, 8, or 9, wherein: It also includes step nine: applying the suborbital vehicle carrying performance optimization method considering flight state constraints described in steps one to eight to suborbital vehicle carrying performance optimization to improve the performance and design efficiency of the suborbital vehicle; According to the sub-orbital vehicle optimization result obtained in step eight, the sub-orbital vehicle maximizes the carrying performance under the constraints of take-off mass, maximum axial overload, maximum normal overload, terminal velocity of reentry section, range of reentry section, terminal lift-drag ratio of reentry section, maximum flight dynamic pressure, terminal velocity inclination angle of reentry section, maximum stagnation heat flux density, and maximum wing leading edge heat flux density; the sub-orbital vehicle carrying performance optimization includes the fields of single-stage horizontal take-off and landing sub-orbital vehicle carrying performance optimization, two-stage horizontal take-off and landing sub-orbital vehicle carrying performance optimization, and vertical take-off and landing sub-orbital vehicle carrying performance optimization.

Images