CN108121856B - Dynamic stability analysis method for full-flight-domain aircraft - Google Patents

Dynamic stability analysis method for full-flight-domain aircraft Download PDF

Info

Publication number
CN108121856B
CN108121856B CN201711277166.8A CN201711277166A CN108121856B CN 108121856 B CN108121856 B CN 108121856B CN 201711277166 A CN201711277166 A CN 201711277166A CN 108121856 B CN108121856 B CN 108121856B
Authority
CN
China
Prior art keywords
equation
aircraft
model
dynamic stability
sampling
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201711277166.8A
Other languages
Chinese (zh)
Other versions
CN108121856A (en
Inventor
张陈安
刘�文
马兴普
雷麦芳
王晓朋
王发民
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Institute of Mechanics of CAS
Original Assignee
Institute of Mechanics of CAS
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Institute of Mechanics of CAS filed Critical Institute of Mechanics of CAS
Priority to CN201711277166.8A priority Critical patent/CN108121856B/en
Publication of CN108121856A publication Critical patent/CN108121856A/en
Application granted granted Critical
Publication of CN108121856B publication Critical patent/CN108121856B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/15—Vehicle, aircraft or watercraft design
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Geometry (AREA)
  • General Physics & Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Aviation & Aerospace Engineering (AREA)
  • Automation & Control Theory (AREA)
  • Feedback Control In General (AREA)
  • Other Investigation Or Analysis Of Materials By Electrical Means (AREA)

Abstract

The invention provides a method for solving the characteristic matrix characteristic root in the state space coupling of an unsteady aerodynamic agent-reduced order model based on efficient self-adaptive sampling and a rigid body dynamic equation to predict the dynamic stability characteristics of the full flight domain, and the coupling dynamic stability characteristics of an aircraft in the full flight domain can be efficiently and accurately acquired. The method can be used for the calculation and analysis of the multi-channel coupling dynamic stability characteristics of aerospace aircrafts, particularly aircrafts with non-axisymmetric complex layouts. The method can greatly reduce the calculation cost while ensuring the prediction precision, thereby being easy to carry out qualitative and quantitative research on the dynamic stability characteristics of the aircraft in the full flight area and obtaining the design direction of guidance significance for the design of the aircraft.

Description

Dynamic stability analysis method for full-flight-domain aircraft
Technical Field
The invention relates to the field of flight mechanics, in particular to a hypersonic aircraft full-flight-domain coupling dynamic stability characteristic analysis method which can meet the precision requirement at low calculation cost.
Background
In the aerospace engineering, for an axisymmetric blunt body layout aircraft, the stability of the aircraft is judged by using a static derivative and a dynamic derivative in China, and no dynamic coupling characteristic appears in the engineering practice. In order to obtain good lift-drag ratio performance, the pneumatic layout of the hypersonic high lift-drag ratio complex layout has flat and non-axisymmetric characteristics, and a dynamic coupling phenomenon can occur in the flying process, so that the problem of failed test flight of the U.S. HTV-2 hypersonic glide aircraft in 2010 is caused.
At present, a common method for researching the coupling dynamic characteristics of an aircraft is computational fluid dynamics/rigid body dynamics (CFD/RBD) coupling numerical solution, in the solution process, the main calculation cost comes from unsteady CFD calculation, under the hypersonic flight condition, a large computer is required to perform calculation for several months to complete the dynamic stability characteristic evaluation in a design state, in the aircraft design process, the dynamic stability characteristics of the aircraft in the whole flight envelope are required to be obtained, and the CFD/RBD coupling of each design point is impractical. The simplified flight mechanics analysis model is constructed by combining the engineering aerodynamic model and the rigid body dynamics model, the calculated amount can be greatly reduced, but a strong viscous interference phenomenon can occur in the hypersonic flight in the near space, so that the precision of most of the simplified aerodynamic models based on the inviscid hypothesis is poor.
In addition, because flight state parameters such as the flight speed, the flight altitude, the flight attitude and the like of the aircraft change in a wide range, when the dynamic characteristic of the whole flight domain of the aircraft is evaluated, only part of state points can be generally calculated or experimented, and the state points which are not involved are generally obtained by interpolation. In this process, the process of selecting the state points for calculation or experiment is called sampling. To ensure interpolation accuracy over the full flight range, it is often necessary to perform intensive sampling over a large range of flight state parameters, and when more flight parameters are involved, it is inevitable to encounter the problem of "dimension cursing", i.e. the number of samples increases geometrically. In addition, the problem of losing important sampling domains randomly is also easy to occur, namely for local areas with violent change of flight mechanical characteristics, the problem that the local interpolation error is too high due to ineffective coverage of random sampling points in the sampling process is solved.
The self-adaptive sampling method can automatically optimize the distribution of sampling points in a large-range and multi-parameter flight domain through the parameter sensitivity of aerodynamic characteristic change, automatically increase the sampling points in a proper region according to the accuracy performance of an evaluation interpolation model, obtain ideal sampling distribution after number theory optimization, and obtain better interpolation accuracy with lower calculation cost. However, in the adaptive sampling process, a large amount of calculation is required to be consumed for estimating the model accuracy, the CFD calculation itself is very expensive for the dynamic stability characteristic research of flight mechanics, the conventional method for directly and randomly extracting the sample points and the CFD results to compare and estimate errors is difficult to perform the global accuracy estimation of the adaptive model at a reasonable calculation cost, and the CFD calculation, especially the CFD calculation related to the unsteady problem, is not high in robustness, so that the usability and efficiency of the adaptive sampling method are not ideal.
Disclosure of Invention
The invention aims to provide an analysis method capable of efficiently and accurately acquiring the coupling dynamic stability characteristics of the full flight domain of an aircraft.
Particularly, the invention provides a dynamic stability analysis method of a full flight domain aircraft, which comprises the following steps:
step 100, determining target precision requirements and calculation upper limits according to aerodynamic characteristics of an aircraft to be researched, and selecting a preset number of initial sample points in a full flight domain range;
at step 110, at the initial sample point x(i)For a given set of motion training signal sequences u (x)(i)K) obtaining a time series y (x) of unsteady aerodynamic and moment coefficients of the forced movement of the aircraft using CFD calculations(i)K), then with u (x)(i)K) and y (x)(i)K) is sample data, and a Kriging agent model is constructed to serve as an initial target agent model;
step 120, constructing a reference agent model, and evaluating the initial sample point to determine the current sampling precision;
step 130, generating new candidate sampling points by using a test design method; calculating and evaluating all candidate sampling points through the target agent model and the reference agent model respectively to obtain precision variation values generated after each candidate sampling point is added into the target agent model and the reference agent model;
140, selecting candidate sampling points meeting the current sampling precision requirement in the precision variation value as next sample points added into the target agent model and the reference agent model for training through a self-adaptive sampling criterion;
150, repeating the steps 130 to 140, then judging whether the precision of the current target agent model reaches the standard and the calculated amount reaches the upper limit or not when repeating every time, and exiting the current cycle when one of the conditions is met, completing the sampling of the current target agent model and determining the training sample of the full-flight-domain unsteady aerodynamic agent-reduced model;
200, generating a large number of design points which are randomly distributed and fill the full flight domain in the full flight domain by using a random sampling method;
step 210, for each design point, obtaining input and output signals of each design point by using Kriging interpolation in a training sample, and determining an unsteady aerodynamic step-down model by using a discrete space formed by each input and output signal;
step 220, converting the unsteady aerodynamic reduced-order model into a state space equation of a continuous space; converting the rigid body power equation into a state equation under a continuous space; performing feedback connection on a state space equation of the unsteady aerodynamic reduced-order model and a state equation of a rigid body dynamic equation to obtain a coupling dynamic stability analysis equation of the current aircraft;
step 230, solving a characteristic matrix characteristic root of the coupled dynamic stability analysis equation, wherein a real part of the characteristic root represents system damping and an imaginary part represents system frequency; when all the real parts of all the characteristic roots are negative, the aircraft at the design point is stable in motion; when the characteristic root of the positive real part appears, the aircraft at the design point is unstable;
and 240, after the dynamic stability characteristics of the aircraft at each design point are obtained through the steps, the dynamic stability characteristics of the aircraft in the whole flight domain can be obtained.
In one embodiment of the present invention, in the step 100, the initial sample points are generated by using O L HS random method, and need to completely cover the whole of the full flight domain space and its boundary vertices.
In one embodiment of the present invention, in the step 140, when the adaptive sampling criterion selects the candidate sampling points, the candidate sampling points whose distances from the existing sampling points are within a specified range and outside the specified range need to be removed.
In one embodiment of the present invention, the functional metric of the adaptive sampling criterion is:
C(ξi)＝(ξi)*ρ(ξi) 1
wherein (ξ)i) Characterizing the reference proxy model and the target proxy model at ξiDeviation of the prediction of (a);
(ξi) The calculation method of (c) is as follows:
Figure BDA0001496798210000031
wherein x(j)Is an existing sample point of the sample, and,
Figure BDA0001496798210000032
and
Figure BDA0001496798210000033
prediction of the reference proxy model and the target proxy model, respectively, when the candidate sample point ξiWith some existing sample point x(j)When in superposition, the zero is forced to eliminate the influence of rounding errors;
ρ(ξi) The expression of (a) is as follows:
Figure BDA0001496798210000034
wherein the molecule is the candidate sampling point ξiThe minimum normalized Euclidean distance between the current sample point and the current sample point, and the denominator is the candidate sample point ξiThe expected closest distance to all existing sample points.
In one embodiment of the invention, the desired nearest neighbor distance is given by:
Figure BDA0001496798210000041
wherein Ns is the number of all existing sample points in the full flight domain, and a is a geometric characterization of the parameter range standardization under study, calculated as follows:
Figure BDA0001496798210000042
x(upper)and x(lower)Upper and lower bounds of the variable dimensions of the parameter space are respectively represented:
Figure BDA0001496798210000043
Figure BDA0001496798210000044
s is Nvar×NvarDimension diagonal matrix, where the ith diagonal element corresponds to the ith variable dimension xlStandard deviation of upper and lower bounds:
Figure BDA0001496798210000045
in the formula oflIs xlMean of upper and lower bounds:
Figure BDA0001496798210000046
in an embodiment of the present invention, the structure of the unsteady aerodynamic reduced-order model in step 210 is:
Figure BDA0001496798210000047
wherein y (k) represents the generalized aerodynamic vector at the kth time step, and u (k-i) represents the generalized displacement input vector at the kth-i time step; a. theiAnd BiIs a coefficient matrix to be identified; na and nb are the delay orders of the output and input, respectively, characterizing the unsteady effects of aerodynamic forces.
In one embodiment of the present invention, when na is 0 and nb is 1, the unsteady aerodynamic reduced-order model structure is essentially a steady model; when na is 0 and nb is more than 1, the unsteady aerodynamic reduced-order model structure is a quasi-steady model in nature; when na is more than 0 and nb is more than 1, the unsteady aerodynamic reduced-order model structure is an unsteady model in nature.
In one embodiment of the present invention, the unsteady aerodynamic reduced order model structure is that the response of the kth time step in the equation 10 is regarded as a linear combination of the first na outputs and nb inputs nearest to the kth time step.
In one embodiment of the present invention, in step 210, the process of converting the unsteady aerodynamic reduced-order model into the state space equation is as follows:
let the modal displacement ξ be u, the modal aerodynamic coefficient faY, define the following state vector:
xa(k)＝[fa(k-1),…,fa(k-na),ξ(k-1),…,ξ(k-nb+1)]T11
equation 10 can be written as:
Figure BDA0001496798210000051
in the formula:
Figure BDA0001496798210000052
Figure BDA0001496798210000053
Figure BDA0001496798210000054
Figure BDA0001496798210000055
the equation 12 is then converted to a continuous system state space equation of equation 13:
Figure BDA0001496798210000056
in one embodiment of the invention, the rigid body equation of power is of the generalized form:
Figure BDA0001496798210000057
wherein M is a generalized mass matrix, Q is a generalized aerodynamic matrix, and ξ is a generalized modal displacement matrix;
defining state variables conditioned on two-degree-of-freedom coupling
Figure BDA0001496798210000058
The rigid body dynamics equation shown in equation 14 can be changed to:
Figure BDA0001496798210000059
in the formula (I), the compound is shown in the specification,
Figure BDA00014967982100000510
Figure BDA00014967982100000511
q represents a free incoming flow pressure;
by applying a discrete space-to-continuous space conversion method, equation 15 can be directly converted into a continuous system state space form of equation 16:
Figure BDA0001496798210000061
in one embodiment of the invention, a state space equation 13 formula of the unsteady aerodynamic reduced order model in the continuous system and a state equation 16 formula of the rigid body power equation in the continuous space are connected in a feedback manner to obtain a coupling dynamic stability analysis equation of the current aircraft, wherein the equation is as follows:
Figure BDA0001496798210000062
aiming at the problem of analyzing the full-flight-domain multi-channel coupling dynamic stability of the aircraft, the invention provides a method for solving a characteristic matrix characteristic root in a state space by coupling an unsteady aerodynamic reduced order model based on a high-efficiency self-adaptive sampling method and a rigid body dynamic equation to predict the full-flight-domain dynamic stability characteristic of the aircraft, which can greatly reduce the calculation cost while ensuring the prediction precision, thereby being easy to carry out qualitative and quantitative researches on the full-flight-domain dynamic stability characteristic of the aircraft and obtaining the design direction with guiding significance on the design of the aircraft.
Drawings
FIG. 1 is a flow chart of the steps of one embodiment of the present invention;
FIG. 2 is a schematic analysis flow diagram of one embodiment of the present invention;
FIG. 3 is a displacement signal diagram of one embodiment of the present invention;
FIG. 4 is an exemplary waverider for one embodiment of the present invention;
FIG. 5 is an adaptive sample point distribution according to an embodiment of the present invention;
FIG. 6 is a full flight domain lateral heading coupled dynamic stability mode profile of an embodiment of the invention.
Detailed Description
As shown in fig. 1 and 2, the method for analyzing the dynamic stability of the full-flight-domain aircraft disclosed by the invention firstly needs to establish a full-flight-domain unsteady aerodynamic agent-reduced order model, then carries out state space coupling dynamic stability analysis aiming at the full-flight domain, and then can obtain the dynamic stability characteristics of the aircraft in the whole flight domain. The process generally comprises the following steps:
step 100, determining target precision requirements and calculation upper limits according to aerodynamic characteristics of an aircraft to be researched, and selecting a preset number of initial sample points in a full flight domain range;
the target precision requirement refers to the precision which is finally realized after the model is subjected to interpolation processing, and the upper limit of the calculated amount is the number of times of adding the candidate sampling points.
Here, theThe initial sample points of (2) should generally contain the boundary vertices of the parameter space (full flight domain range) to avoid extrapolation in proxy model prediction. But for very high parameter dimensions, all vertices are counted to the initial sample point, also encountering the "cursing" of the dimension. For example, for a 12-dimensional parameter space, the total number of vertices is 4096 (i.e., 2)12) In this case, a random method such as Optimized-L where Hypercube Sampling (O L HS) is usually used to generate a set of sample points that can substantially completely cover the whole parameter space for sample initialization.
At step 110, at the initial sample point x(i)For a given set of motion training signal sequences u (x)(i)K) obtaining a time series y (x) of unsteady aerodynamic and moment coefficients of the forced movement of the aircraft using CFD calculations(i)K), then with u (x)(i)K) and y (x)(i)K) is sample data, and a Kriging agent model is constructed to serve as an initial target agent model;
the motion training signal refers to a rigid motion time-varying process that an aircraft meets a certain frequency characteristic, and usually, a '3211' signal can meet the frequency characteristic requirement. For example, to predict the roll/yaw coupled dynamic stability characteristics of the aircraft, the aerodynamic forces and moments of the aircraft during simultaneous roll/yaw coupled motions are solved through CFD. The input signal is the '3211' signal shown in fig. 2, wherein 'x 1' represents the roll angular displacement, 'x 2' represents the yaw angular displacement, and the calculated aerodynamic force and moment are the output.
Here, although an aircraft is used as the target, an airfoil may be used as the target.
Step 120, constructing a reference agent model, and evaluating the initial sample point to determine the current sampling precision;
in order to reduce the computation amount of the unsteady CFD, in addition to constructing a target proxy Model (DM), a Reference proxy Model (RM) is also constructed by other interpolation methods for assisting sampling. In each round of adaptive sampling, assuming that the prediction of the RM at the sampling candidate point is true, the prediction of the RM can be used to estimate the prediction error of the DM at the same point. In the whole sampling process, the unsteady CFD calculation is only used for training new sample points, and the test sample points are not required to be calculated, so that the calculation consumption can be greatly saved, and the modeling efficiency of the target agent model is improved.
Step 130, generating new candidate sampling points by using a test design method; calculating and evaluating all candidate sampling points through the target agent model and the reference agent model respectively to obtain precision variation values generated after each candidate sampling point is added into the target agent model and the reference agent model;
the existing (Design of experience) DoE method is adopted to randomly generate a predetermined number of candidate sampling points, and the number of specific candidate sampling points can be set randomly according to needs.
In the step, each candidate sampling point is respectively substituted into the DM for calculation so as to obtain an influence value on the model precision after the corresponding candidate sampling point is added into the model in advance, then the RM evaluates the influence value of the candidate sampling point, and finally the sequence of all the candidate sampling points taking the influence value as the arrangement sequence can be obtained.
When proxy model verification is performed, there are various error metric indexes, such as Nash-Sutcliffe efficiency metric (NSEM), which is defined as:
Figure BDA0001496798210000081
in the formula (I), the compound is shown in the specification,
Figure BDA0001496798210000082
indicates the predicted value, f (x)i) A true value is shown in the table of values,
Figure BDA0001496798210000083
denotes f (x)i) And (4) average value.
140, selecting candidate sampling points meeting the current sampling precision requirement in the precision variation value as next sample points added into the target agent model and the reference agent model for training through a self-adaptive sampling criterion;
the adaptive Sampling criterion (Sampling criterion) in this step has an important influence on the sample point distribution and the prediction accuracy of the proxy model. In sampling optimization, the current view is unified, so that the sample points can not only fill the space, but also meet the requirement of local reinforcement. The full space is the measurement of the spatial distribution of the sample points, and the full parameter space of the sample points is beneficial to improving the global prediction precision of the proxy model; and local reinforcement is a measurement on function response, and a sample point is locally reinforced in some important sampling domains, so that the capability of the proxy model for predicting the nonlinear characteristics of the function response can be enhanced.
When the self-adaptive sampling criterion selects the candidate sampling points, the candidate sampling points with the distance between the existing sampling points within a specified range and the distance between the existing sampling points and the candidate sampling points outside the specified range need to be removed. For the adaptive selection of the sample points, the key point is to adjust the geometric constraint of the distribution of the sample points, not only to avoid that a new candidate sampling point is too close to an existing sampling point, but also to prevent the new candidate sampling point from being too far away, when the sample points are gathered in a certain or a plurality of local sampling domains, the overfitting of the proxy model can be caused, and when the sample points are too far away, the capability of the proxy model for capturing the nonlinear characteristic of the function response can be weakened.
In addition, when the candidate sampling point is overlapped with the existing sampling point, the candidate sampling point is rejected.
Therefore, in this step, the candidate sampling point selected by the adaptive sampling criterion does not necessarily have the greatest influence on the accuracy of the current model, but is distributed most reasonably.
150, repeating the steps 130 to 140, then judging whether the precision of the current target agent model reaches the standard and the calculated amount reaches the upper limit or not when repeating every time, and exiting the current cycle when one of the conditions is met, completing the sampling of the current target agent model and determining the training sample of the full-flight-domain unsteady aerodynamic agent-reduced model;
when the sampling points are added, the steps are repeated, and the optimal sampling point in the next batch of candidate sampling points can be selected according to the self-adaptive sampling criterion on the basis of adding new sampling points each time.
When the precision of the added candidate sampling point meets the set precision requirement, the supplement of the sampling point can be stopped. Or under the condition that the precision does not meet the requirement, but the number of the increased sampling points meets the requirement, the supplement of the sampling points can be stopped, and the setting can avoid using excessive calculation time.
After the training sample is established, for any new design point, the input and output signals of the design point can be obtained through Kriging interpolation.
In the foregoing steps, the function metric of the specific adaptive sampling criterion is:
C(ξi)＝(ξi)*ρ(ξi) (1)
wherein (ξ)i) Characterization RM and DM at ξiThe deviation of the prediction at (d) is called the Deviation Function (DF).
The method for updating the candidate sample point is as follows:
1) evaluating the value of the adaptive sampling criterion C at all Ncand pre-generated candidate sampling points;
2) and the candidate sampling points are arranged in a descending order according to the value of C obtained in the previous step, wherein the first candidate sampling point is a state point observing the maximum C value, and is selected as a new sample point to wait for the training of the original model.
3) In each round of adaptive sampling, the above steps are repeated until there are number of additional points per circle generated.
(ξi) The calculation method of (c) is as follows:
Figure BDA0001496798210000091
wherein x(j)Is an existing sample point of the sample, and,
Figure BDA0001496798210000092
and
Figure BDA0001496798210000093
prediction of RM and DM, respectively, as candidate sample point ξiWith some existing sample point x(j)When in superposition, the zero is forced to eliminate the influence of rounding errors;
ρ(ξi) The expression of (a) is as follows:
Figure BDA0001496798210000094
wherein the molecule is the candidate sampling point ξiThe minimum normalized Euclidean distance between the current sample point and the current sample point, and the denominator is the candidate sample point ξiThe expected closest distance to all existing sample points.
ρ(ξi) The ratio of the minimum distance to the expected distance is a measure of the spatial distribution pattern of the candidate sampling points and the existing sampling points, and the non-sampling region and the sampling region can be effectively distinguished. After the parameter range is determined, the desired distance rexpOnly with respect to the total number of sample points Ns. Expected distance r in adaptive sampling processexpThe ratio ρ (ξ) of the minimum distance to the desired distance is also updated continuously as the total number of sample points Ns increasesi) For adaptive sampling of sample points, the significance of this parameter is to adjust the geometric constraints of the distribution of sample points, a large q means that the distance between sample points is limited to be increased, so that the sample points can be effectively prevented from being clustered in some local sampling domain or domains, resulting in overfitting of the proxy modeli) The separation/clustering characteristics of the spatial distribution of sample points are known and are therefore called separation functions.
Further, the desired closest proximity distance in the preceding step may be given by:
Figure BDA0001496798210000101
wherein Ns is the number of all existing sample points in the parameter space, and a is a geometric characterization of the parameter range standardization under study, calculated as follows:
Figure BDA0001496798210000102
x(upper)and x(lower)Upper and lower bounds of the variable dimensions of the parameter space are respectively represented:
Figure BDA0001496798210000103
Figure BDA0001496798210000104
s is Nvar×NvarDimension diagonal matrix, where the ith diagonal element corresponds to the ith variable dimension xlStandard deviation of upper and lower bounds:
Figure BDA0001496798210000105
in the formula oflIs xlMean of upper and lower bounds:
Figure BDA0001496798210000106
200, generating a large number of design points which are randomly distributed and fill the full flight domain in the full flight domain by using a random sampling method;
step 210, for each design point, obtaining input and output signals of each design point by using Kriging interpolation in a training sample, and determining an unsteady aerodynamic step-down model by using a discrete space formed by each input and output signal;
Figure BDA0001496798210000107
wherein y (k) represents the generalized pneumatic force vector at the k time step, and u (k-i) represents the generalized displacement input vector at the k-i time step; a. theiAnd BiIs a coefficient matrix to be identified; na and nb are the delay orders of the output and input, respectively, characterizing the unsteady effects of aerodynamic forces.
Sampling least square identification solution undetermined coefficient matrix AiAnd BiThe required unsteady aerodynamic force model is obtained, and according to experience, the accuracy requirement can be met by taking na and nb as 3-5 generally.
Step 220, converting the unsteady aerodynamic reduced-order model into a state space equation of a continuous space; converting the rigid body power equation into a state equation under a continuous space; performing feedback connection on a state space equation of the unsteady aerodynamic reduced-order model and a state equation of a rigid body dynamic equation to obtain a coupling dynamic stability analysis equation of the current aircraft; the transformation process of the unsteady aerodynamic reduced-order model is as follows:
let the modal displacement ξ be u, the modal aerodynamic coefficient faY, define the following state vector:
xa(k)＝[fa(k-1),…,fa(k-na),ξ(k-1),…,ξ(k-nb+1)]T(11)
equation (10) can be written as:
Figure BDA0001496798210000111
in the formula:
Figure BDA0001496798210000112
Figure BDA0001496798210000113
Figure BDA0001496798210000114
Figure BDA0001496798210000115
fa0the static aerodynamic force coefficient does not influence the dynamic stability of the aircraft, and can be ignored in the dynamic stability analysis;
the method for converting the discrete space into the continuous space can directly convert (12) into a continuous system state space form of the formula (13):
Figure BDA0001496798210000116
the transformation process of the rigid body power equation transformation is as follows:
the generalized form of the rigid body equation is:
Figure BDA0001496798210000117
wherein M is a generalized mass matrix, Q is a generalized aerodynamic matrix, and ξ is a generalized modal displacement matrix;
defining state variables with two degree-of-freedom coupling as an example (e.g. roll/yaw coupling)
Figure BDA0001496798210000118
The rigid body dynamics system represented by equation (14) may be changed to:
Figure BDA0001496798210000121
in the formula (I), the compound is shown in the specification,
Figure BDA0001496798210000122
Figure BDA0001496798210000123
q represents a free incoming flow pressure.
The method for converting the discrete space into the continuous space can be used for directly converting (15) into a continuous system state space form of formula (16):
Figure BDA0001496798210000124
the current coupled dynamic stability analysis equation of the aircraft is as follows:
Figure BDA0001496798210000125
step 230, solving a characteristic matrix characteristic root of the coupled dynamic stability analysis equation, wherein a real part of the characteristic root represents system damping and an imaginary part represents system frequency; when all the real parts of all the characteristic roots are negative, the aircraft at the design point is stable in motion; when the characteristic root of the positive real part appears, the aircraft at the design point is unstable;
and 240, after the dynamic stability characteristics of the aircraft at each design point are obtained through the steps, the dynamic stability characteristics of the aircraft in the whole flight domain can be obtained.
Exemplary embodiments
The method for analyzing the coupling dynamic stability characteristics of the aircraft is utilized to research the lateral direction coupling dynamic stability characteristics of the waverider aircraft in the full flight domain shown in figure 4. flight design parameters are selected from Mach number M, altitude H and attack angle α, the change range of the Mach number is 5-20, the change range of the altitude is 25-70 km, and the change range of the attack angle α is-10-20 degrees.
First, the number of initialization sample points of the proxy-reduced model is N0Other parameters of the adaptive sampling are q 2 and N30s1＝20、N m10, p 1 and Ncand6000. Model output dimension N in this examplerep3, the 3 sample points are updated in each sampling iteration for each response, so the total number of sample points for each update is Nuppc9. To save the time consuming non-stationary CFD calculations, sampling is terminated when the total number of sample points exceeds 200 or the NSEM (Nash-Sutcliffe efficiency metric) estimates for both responses in the adaptive sampling reach or exceed 0.9. The final training sample point number of this example is N s30+9 × 16, 174, where the number of sample points updated by the adaptive method is 144,the sample point distribution is shown in fig. 5.
Then, the '3211' displacement signal shown in fig. 3 is used as an input signal, and each training sample point is substituted into the CFD solution to obtain an aerodynamic output corresponding to the input signal, including a roll moment coefficient and a yaw moment coefficient.
Then, 11880 design sample points are selected in the flight domain design space through Monte Carlo sampling, the input and output signals of each design point are obtained through a proxy model and are substituted into the formula (10), and the coefficient matrix A is obtained through least square identificationiAnd BiAnd determining the unsteady aerodynamic order-reducing model.
Finally, through the implementation of steps 210 to 260, the characteristic matrix characteristic root of each design point can be obtained, and the lateral heading coupled dynamic stability characteristic of the aircraft in the whole flight domain is also obtained, the distribution of the lateral heading coupled dynamic stability mode in the whole flight domain is shown in fig. 6, wherein the reference attack angle of the graph (a) is α00 DEG, the reference angle of attack in the diagram (b) is α0In the figure, "Dutch Roll Divergence" indicates the Dutch Roll Divergence mode, "Dutch Roll Divergence" indicates the Dutch Roll Convergence mode, and "Static Divergence" indicates the Roll Divergence mode, 10 °.
The above embodiment shows the whole process of analyzing the coupling dynamic stability characteristics of the aircraft in the full flight domain by the full flight domain aircraft dynamic stability analysis method provided by the patent.
Thus, it should be understood by those skilled in the art that while exemplary embodiments of the present invention have been illustrated and described in detail herein, many other variations or modifications which are consistent with the principles of the invention may be directly determined or derived from the disclosure of the present invention without departing from the spirit and scope of the invention. Accordingly, the scope of the invention should be understood and interpreted to cover all such other variations or modifications.

Claims (10)

1. A full flight domain aircraft dynamic stability analysis method is characterized by comprising the following steps:
step 100, determining target precision requirements and calculation upper limits according to aerodynamic characteristics of an aircraft to be researched, and selecting a preset number of initial sample points in a full flight domain range;
at step 110, at the initial sample point x(i)For a given set of motion training signal sequences u (x)(i)K) obtaining a time series y (x) of unsteady aerodynamic and moment coefficients of the forced movement of the aircraft using CFD calculations(i)K), then with u (x)(i)K) and y (x)(i)K) is sample data, and a Kriging agent model is constructed to serve as an initial target agent model;
step 120, constructing a reference agent model, and evaluating the initial sample point to determine the current sampling precision;
step 130, generating new candidate sampling points by using a test design method; calculating and evaluating all candidate sampling points through the target agent model and the reference agent model respectively to obtain precision variation values generated after each candidate sampling point is added into the target agent model and the reference agent model;
140, selecting candidate sampling points meeting the current sampling precision requirement in the precision variation value as next sample points added into the target agent model and the reference agent model for training through a self-adaptive sampling criterion;
150, repeating the steps 130 to 140, then judging whether the precision of the current target agent model reaches the standard and the calculated amount reaches the upper limit or not when repeating every time, and exiting the current cycle when one of the conditions is met, completing the sampling of the current target agent model and determining the training sample of the full-flight-domain unsteady aerodynamic agent-reduced model;
200, generating a large number of design points which are randomly distributed and fill the full flight domain in the full flight domain by using a random sampling method;
step 210, for each design point, obtaining input and output signals of each design point by using Kriging interpolation in a training sample, and determining an unsteady aerodynamic step-down model by using a discrete space formed by each input and output signal;
step 220, converting the unsteady aerodynamic reduced-order model into a state space equation of a continuous space; converting the rigid body power equation into a state equation under a continuous space; performing feedback connection on a state space equation of the unsteady aerodynamic reduced-order model and a state equation of a rigid body dynamic equation to obtain a coupling dynamic stability analysis equation of the current aircraft;
step 230, solving a characteristic matrix characteristic root of the coupled dynamic stability analysis equation, wherein a real part of the characteristic root represents system damping and an imaginary part represents system frequency; when all the real parts of all the characteristic roots are negative, the aircraft at the design point is stable in motion; when the characteristic root of the positive real part appears, the aircraft at the design point is unstable;
and 240, after the dynamic stability characteristics of the aircraft at each design point are obtained through the steps 210 to 230, the dynamic stability characteristics of the aircraft in the whole flight domain can be obtained.
2. The full flight domain aircraft dynamic stability analysis method of claim 1,
in the step 100, initial sample points are generated by using an O L HS random method, and need to completely cover the whole flight domain space and its boundary vertices.
3. The full flight domain aircraft dynamic stability analysis method of claim 1,
in step 140, when the adaptive sampling criterion selects the candidate sampling points, the candidate sampling points whose distances from the existing sampling points are within the specified range and outside the specified range need to be removed.
4. The full flight domain aircraft dynamic stability analysis method of claim 3,
the functional metric of the adaptive sampling criterion is:
C(ξi)＝(ξi)*ρ(ξi) 1
wherein (ξ)i) Characterizing the reference proxy model and the target proxy model at ξiDeviation of the prediction of (a);
(ξi) The calculation method of (c) is as follows:
Figure FDA0002448129230000021
wherein x(j)Is an existing sample point of the sample, and,
Figure FDA0002448129230000022
and
Figure FDA0002448129230000023
prediction of the reference proxy model and the target proxy model, respectively, when the candidate sample point ξiWith some existing sample point x(j)When in superposition, the zero is forced to eliminate the influence of rounding errors;
ρ(ξi) The expression of (a) is as follows:
Figure FDA0002448129230000031
wherein the molecule is the candidate sampling point ξiThe minimum normalized Euclidean distance between the current sample point and the current sample point, and the denominator is the candidate sample point ξiThe expected closest distance to all existing sample points.
5. The full flight domain aircraft dynamic stability analysis method of claim 4,
the desired nearest neighbor distance is given by:
Figure FDA0002448129230000032
wherein Ns is the number of all existing sample points in the full flight domain, and a is a geometric characterization of the parameter range standardization under study, calculated as follows:
Figure FDA0002448129230000033
x(upper)and x(lower)Upper and lower bounds of the variable dimensions of the parameter space are respectively represented:
Figure FDA0002448129230000034
Figure FDA0002448129230000035
s is Nvar×NvarDimension diagonal matrix, where the ith diagonal element corresponds to the ith variable dimension xlStandard deviation of upper and lower bounds:
Figure FDA0002448129230000036
in the formula oflIs xlMean of upper and lower bounds:
Figure FDA0002448129230000037
6. the full flight domain aircraft dynamic stability analysis method of claim 1,
the structure of the unsteady aerodynamic reduced-order model in the step 210 is as follows:
Figure FDA0002448129230000038
wherein y (k) represents the generalized aerodynamic vector at the kth time step, and u (k-i) represents the generalized displacement input vector at the kth-i time step; a. theiAnd BiIs a coefficient matrix to be identified; na and nb are the delay orders of the output and input, respectively, characterizing the unsteady effects of aerodynamic forces.
7. The full flight domain aircraft dynamic stability analysis method of claim 6,
when na is 0 and nb is 1, the unsteady aerodynamic reduced-order model structure is essentially a steady model; when na is 0 and nb is more than 1, the unsteady aerodynamic reduced-order model structure is a quasi-steady model in nature; when na is more than 0 and nb is more than 1, the unsteady aerodynamic reduced-order model structure is an unsteady model in nature.
8. The full flight domain aircraft dynamic stability analysis method of claim 7,
in step 210, the process of converting the unsteady aerodynamic reduced-order model into the state space equation is as follows:
let the modal displacement ξ be u, the modal aerodynamic coefficient faY, define the following state vector:
xa(k)＝[fa(k-1),…,fa(k-na),ξ(k-1),…,ξ(k-nb+1)]T11
then equation 10 can be written as:
Figure FDA0002448129230000041
in the formula:
Figure FDA0002448129230000042
Figure FDA0002448129230000043
Figure FDA0002448129230000044
Figure FDA0002448129230000045
equation 12 is then converted to equation 13 for the state space equation for a continuous system:
Figure FDA0002448129230000051
9. the full flight domain aircraft dynamic stability analysis method of claim 8,
the generalized form of the rigid body equation is:
Figure FDA0002448129230000052
wherein M is a generalized mass matrix, Q is a generalized aerodynamic matrix, and ξ is a generalized modal displacement matrix;
defining state variables conditioned on two-degree-of-freedom coupling
Figure FDA0002448129230000053
The rigid body dynamics equation shown in equation 14 can be changed to:
Figure FDA0002448129230000054
in the formula (I), the compound is shown in the specification,
Figure FDA0002448129230000055
Figure FDA0002448129230000056
q represents a free incoming flow pressure;
by applying a discrete space-to-continuous space conversion method, equation 15 can be directly converted into a continuous system state space form of equation 16:
Figure FDA0002448129230000057
10. the method for analyzing the dynamic stability of the full-flight-domain aircraft according to claim 9,
and performing feedback connection on a state space equation 13 formula of the unsteady aerodynamic reduced-order model in a continuous system and a state equation 16 formula of the rigid body dynamic equation in a continuous space to obtain a coupling dynamic stability analysis equation of the current aircraft, wherein the equation is as follows:
Figure FDA0002448129230000058
CN201711277166.8A 2017-12-06 2017-12-06 Dynamic stability analysis method for full-flight-domain aircraft Active CN108121856B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201711277166.8A CN108121856B (en) 2017-12-06 2017-12-06 Dynamic stability analysis method for full-flight-domain aircraft

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201711277166.8A CN108121856B (en) 2017-12-06 2017-12-06 Dynamic stability analysis method for full-flight-domain aircraft

Publications (2)

Publication Number Publication Date
CN108121856A CN108121856A (en) 2018-06-05
CN108121856B true CN108121856B (en) 2020-08-04

Family

ID=62229823

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201711277166.8A Active CN108121856B (en) 2017-12-06 2017-12-06 Dynamic stability analysis method for full-flight-domain aircraft

Country Status (1)

Country Link
CN (1) CN108121856B (en)

Families Citing this family (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109753690A (en) * 2018-12-10 2019-05-14 西北工业大学 Nonlinear unsteady aerodynamics order reducing method based on Fluid Mechanics Computation
CN109634309B (en) * 2019-02-21 2024-03-26 南京晓庄学院 Autonomous obstacle avoidance system and method for aircraft and aircraft
CN111338364B (en) * 2019-11-21 2021-09-21 浙江大学 High-precision controller for optimizing trajectory of hypersonic aerocraft with quick response
CN111159942B (en) * 2019-12-26 2023-09-15 北京电子工程总体研究所 Method for calculating rolling damping moment of winged aircraft based on steady simulation
CN112699624B (en) * 2021-03-24 2021-06-22 南京信息工程大学 Trajectory calculation method under severe meteorological conditions
CN112711815B (en) * 2021-03-25 2021-06-25 中国科学院自动化研究所 Aircraft modeling and model characteristic analysis system
CN113935256B (en) * 2021-09-27 2024-07-12 北京理工大学 High-dimensional complex aircraft system reduced order characterization method based on error correction
CN113609600B (en) * 2021-10-11 2021-12-14 中国空气动力研究与发展中心计算空气动力研究所 Multi-body separation compatibility measurement and characterization method suitable for aircraft
CN114323552B (en) * 2021-11-18 2022-10-21 厦门大学 Method for judging stability of water entering and exiting from cross-medium navigation body

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102682173A (en) * 2012-05-13 2012-09-19 北京理工大学 Optimization design method based on self-adaptive radial basis function surrogate model for aircraft
CN107220403A (en) * 2017-04-20 2017-09-29 南京航空航天大学 The control association modeling method of aircraft Elastic mode

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102682173A (en) * 2012-05-13 2012-09-19 北京理工大学 Optimization design method based on self-adaptive radial basis function surrogate model for aircraft
CN107220403A (en) * 2017-04-20 2017-09-29 南京航空航天大学 The control association modeling method of aircraft Elastic mode

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
On the generation of flight dynamics aerodynamic tables by computational fluid dynamics;A. Da Ronch et al.;《Progress in Aerospace Sciences》;20111013;全文 *
基于 Kriging 代理模型的飞行器结构刚度气动优化设计;聂雪媛 等;《气体物理》;20170331;全文 *
基于计算试验设计与代理模型的飞行器近似优化策略探讨;龙腾 等;《机械工程学报》;20160731;全文 *

Also Published As

Publication number Publication date
CN108121856A (en) 2018-06-05

Similar Documents

Publication Publication Date Title
CN108121856B (en) Dynamic stability analysis method for full-flight-domain aircraft
CN105843073B (en) A kind of wing structure aeroelastic stability analysis method not knowing depression of order based on aerodynamic force
CN112016167B (en) Aircraft aerodynamic shape design method and system based on simulation and optimization coupling
CN105975645B (en) A kind of aircraft flow field of region containing shock wave quick calculation method based on multistep
Silva et al. Development of reduced-order models for aeroelastic analysis and flutter prediction using the CFL3Dv6. 0 code
Crowell et al. Model reduction of computational aerothermodynamics for hypersonic aerothermoelasticity
Timme et al. Transonic aeroelastic stability analysis using a Kriging-based Schur complement formulation
CN104866692A (en) Aircraft multi-objective optimization method based on self-adaptive agent model
CN114065662B (en) Method suitable for rapidly predicting airfoil flow field with variable grid topology
CN112414668B (en) Wind tunnel test data static bomb correction method, device, equipment and medium
CN104794332B (en) A kind of Uncertainty Analysis Method of skyscraper wind-excited responese analysis model
Güner et al. Unsteady aerodynamic modeling methodology based on dynamic mode interpolation for transonic flutter calculations
CN114564787A (en) Bayesian optimization method, device and storage medium for target-related airfoil design
CN106874561B (en) Multidisciplinary uncertainty propagation analysis method based on Newton iteration
EP2924598A1 (en) A method for determining a structural response of a flow body to an atmospheric disturbance
CN116628854A (en) Wing section aerodynamic characteristic prediction method, system, electronic equipment and storage medium
Apponsah et al. Aerodynamic shape optimization for unsteady flows: some benchmark problems
CN115774900A (en) Instruction robustness optimization design method for variable configuration aircraft under uncertain conditions
CN111177855B (en) Pneumatic structure solving method and system in global aeroelasticity optimization
CN105260498B (en) A kind of large size civil aircraft wing variable camber design method
Amrit et al. Efficient multi-objective aerodynamic optimization by design space dimension reduction and co-kriging
Clarich et al. Reliability-based design optimization applying polynomial chaos expansion: theory and Applications
Barrett et al. Airfoil design and optimization using multi-fidelity analysis and embedded inverse design
Shu et al. Parametric Aeroelastic Reduced-Order Modeling with Hyperparameter Optimization for Flutter Analysis
Padovan et al. Multi objective robust design optimization of airfoils in transonic field

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant