CN108121856A - A kind of full flight domain aerocraft dynamic stability analysis method - Google Patents
A kind of full flight domain aerocraft dynamic stability analysis method Download PDFInfo
- Publication number
- CN108121856A CN108121856A CN201711277166.8A CN201711277166A CN108121856A CN 108121856 A CN108121856 A CN 108121856A CN 201711277166 A CN201711277166 A CN 201711277166A CN 108121856 A CN108121856 A CN 108121856A
- Authority
- CN
- China
- Prior art keywords
- mrow
- msub
- mtd
- mtr
- msup
- 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.)
- Granted
Links
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
Abstract
The present invention provides it is a kind of the unsteady aerodynamic force based on efficient self-adapted sampling acted on behalf of into reduced-order model and Dynamical Equations of Rigid Body solve in State-Space Coupling the method for eigenmatrix characteristic root carry out full flight domain dynamic stability feature and predict, coupled dynamic stability feature that can be with the acquisition aircraft of efficiently and accurately in full flight domain.Available for the calculating analysis to aerospace flight vehicle, particularly the multi-channel coupling dynamic stability characteristic of non-axis symmetry complex topology aircraft.Calculating cost can be greatly decreased while precision of prediction is ensured, so as to be easy to carry out the qualitative and quantitative study of aircraft dynamic stability characteristic in full flight domain, obtain the design direction for having directive significance to Flight Vehicle Design.
Description
Technical field
The present invention relates to flight mechanics fields, can meet required precision more particularly to a kind of with relatively low calculating cost
The full flight domain coupled dynamic stability characteristic analysis method of hypersonic aircraft.
Background technology
In aerospace engineering, for axial symmetry blunt body layout aircraft, China is judged with static derivative and dynamic derivative
There is not dynamical couple characteristic in engineering practice in stability of aircraft.And hypersonic high lift-drag ratio complex topology in order to
Good lift resistance ratio performance is obtained, aerodynamic arrangement is presented flat, non-axis symmetry feature, dynamic coupling is present in flight course
Phenomenon is closed, this problem results in taking a flight test unsuccessfully for the hypersonic glide vehicles of U.S. HTV-2 in 2010.
At present, a kind of common method of exploratory flight device coupling dynamic characteristic is Fluid Mechanics Computation/dynamics of rigid bodies
(CFD/RBD) Coupled Numerical solves, and in solution procedure, main calculating, which expends, comes from unsteady CFD calculating, in high ultrasound
, it is necessary to which the calculating of some months is carried out using mainframe computer could complete dynamic stability under a design point under fast flying condition
Characteristic evaluation, during Flight Vehicle Design, it is also necessary to obtain its dynamic stability characteristic in entire flight envelope, be set to each
It is unpractical that enumeration, which carries out CFD/RBD couplings,.Pass through incorporation engineering Aerodynamic Model and Rigid Body Dynamics Model structure letter
The flight mechanics analysis model of change can substantially reduce calculation amount, but be present near space hypersonic flight strong
Sticky interference phenomenon so that most at present based on poor without the viscous simplification Aerodynamic Model precision assumed.
Further, since the flight status parameters such as vehicle flight speeds, flying height, flight attitude are in wide variation,
In the full flight domain dynamic characteristic assessment of progress aircraft, generally partial status point therein can only be calculated or tested,
The state point being not involved with then generally is obtained by interpolation.In this process, the state point for being calculated or being tested is chosen
Process be referred to as sample.In order to ensure the interpolation precision in full flight domain, it usually needs in large-scale flight status parameter
Intensive sampling is carried out, when being related to more flight parameter, will inevitably meet with " dimension curse " problem, i.e. number of samples is in
Geometric growth.The problem of being additionally susceptible to random loss important sampling domain, i.e., it is acute for flight mechanics characteristic variations
Strong regional area, the stochastical sampling point covering not yet in effect in sampling process and the problem of cause local interpolation error excessively high.
The sensitivity to parameter that adaptively sampled method can be changed by aerodynamic characteristic is realized to fly a wide range of, multi-parameter
Row automatically optimizes sampling point distributions in domain, is increased automatically in appropriate area according to assessment interpolation model precision performance and sampled
Point obtains preferable sample distribution after number theory optimizes, and preferable interpolation precision is obtained with relatively low calculating cost.So
And during adaptively sampled, assessment models precision needs to expend big calculation amount, and flight mechanics dynamic stability characteristic is ground
For studying carefully, cost of CFD calculating itself is very high, and carrying out comparison with CFD results with traditional directly randomly drawing sample point comments
Estimate the method for error, it is difficult to adaptive model global precision assessment is carried out in the case where reasonably calculating cost, and since CFD is counted
It calculates, more particularly to the CFD calculating of Non―steady Problem, robustness itself is not high, causes the availability and effect of adaptively sampled method
Rate is undesirable.
The content of the invention
The invention aims to provide one kind can efficiently and accurately obtain the full flight domain coupled dynamic stability of aircraft
The analysis method of feature.
Particularly, the present invention provides a kind of full flight domain aerocraft dynamic stability analysis method, includes the following steps：
Step 100, for the Aerodynamic Characteristics of required research, determine in aimed at precision requirement and calculation amount
Limit, and choose the initial sample point of predetermined quantity in full flight domain scope；
Step 110, in initial sample point x(i)Place, the training signal sequence u (x given to one group(i), k) and it uses
CFD calculates the when ordinal series y (x of the unsteady aerodynamic force for obtaining aircraft forced movement and torque coefficient(i), k), then with u
(x(i), k) and y (x(i), k) and it is sample data, structure Kriging agent models are as initial target proxy model；
Step 120, structure is assessed initial sample point, with reference to agent model to determine present sample precision；
Step 130, new candidate's sampled point is generated using test design method；And pass through target proxy model and with reference to generation
Reason model is respectively calculated and assesses to all candidate's sampled points, and target proxy model is added in obtain each candidate's sampled point
With the precision changing value generated after reference agent model；
Step 140, the candidate for meeting present sample required precision in precision changing value is chosen by adaptively sampled criterion
Sample point of the sampled point as next addition target proxy model and with reference to agent model is trained；
Step 150, step 130 is repeated to 140, then whether judges current goal agent model precision in each repeat
Whether up to standard and calculation amount reaches the upper limit, and previous cycle is exited when meeting one of condition, completes current goal agency
The sampling of model and the training sample that full flight domain unsteady aerodynamic force agency-reduced-order model is determined with this；
Step 200, a large amount of random distributions, the design full of full flight domain are generated in full flight domain with arbitrary sampling method
Point；
Step 210, to each design point, the input for obtaining each design point using Kriging interpolation in training sample is defeated
Go out signal, the discrete space that each input/output signal is formed is recycled to determine unsteady aerodynamic force reduced-order model；
Step 220, unsteady aerodynamic force reduced-order model is converted into the state space equation of continuous space；By rigid body dynamic
The equations turned state equation under continuous space；By the state space equation and rigid body dynamic of unsteady aerodynamic force reduced-order model
The state equation of equation carries out feedback link to get to the coupled dynamic stability analysis equation of current flight device；
Step 230, the eigenmatrix characteristic root of coupled dynamic stability analysis equation, the characterization system resistance of characteristic root real part are solved
Buddhist nun, imaginary part characterization system frequency；Wherein, when all characteristic root real parts are all negative, the aircraft dynamic stability of the design point；When
When there is the characteristic root of positive real part, the aircraft of the design point moves unstable；
Step 240, after the dynamic stability feature that the aircraft of each design point is obtained by above-mentioned steps, you can flown
Dynamic stability feature of the row device in entire flight domain.
In an embodiment of the invention, in the step 100, initial sample point uses OLHS random devices
Generation, and need the entire full flight domain space of complete covering and its border vertices.
In an embodiment of the invention, in the step 140, the adaptively sampled criterion is choosing candidate
During sampled point, it is necessary to reject relatively have sampled point distance within the specified range with candidate's sampled point outside specified range.
In an embodiment of the invention, the function measurement of the adaptively sampled criterion is：
C(ξi)=δ (ξi)*ρ(ξi) 1
Wherein, δ (ξi) characterize with reference to agent model and target proxy model in ξiThe deviation of the prediction at place；
δ(ξi) calculation it is as follows：
Wherein x(j)It is existing sample point,WithIt is with reference to agent model and target proxy respectively
The prediction of model；As candidate's sampled point ξiHas sample point x with certain(j)During coincidence, δ is forced zero setting, to eliminate rounding error
It influences；
ρ(ξi) expression formula it is as follows：
Wherein, molecule is candidate's sampled point ξiThe minimum standardization Euclidean distance between all existing sample points；Denominator
It is candidate's sampled point ξiThe closest distance of expectation between all existing sample points.
In an embodiment of the invention, it is described it is expected that closest distance is given by：
Wherein Ns is all existing sample point numbers in full flight domain, and A is the geometry of studied parameter area standardization
Characterization, calculates as follows：
x(upper)And x(lower)The upper and lower border of each variable dimension in expression parameter space respectively：
S is Nvar×NvarDiagonal matrix is tieed up, wherein l-th of diagonal entry corresponds to l-th of variable dimension xlIt is upper following
The standard deviation on boundary：
μ in formulalIt is xlThe average of up-and-down boundary：
In an embodiment of the invention, the unsteady aerodynamic force reduced-order model structure in the step 210 is：
Y (k) represents the pneumatic force vector of broad sense of kth time step in formula, and u (k-i) represents the broad sense position of-i time steps of kth
Move input vector；AiAnd BiFor coefficient matrix to be identified；Na and nb is respectively the delay exponent number for exporting and inputting, and characterization is pneumatic
The unsteady aerodynamic effect of power.
In an embodiment of the invention, as na=0 and nb=1, the unsteady aerodynamic force reduced-order model knot
Structure is substantially a kind of steady model；As na=0 and nb ＞ 1, it is on the unsteady aerodynamic force reduced-order model structural nature
A kind of quasi-steady model；It is a kind of non-fixed as na ＞ 0 and nb ＞ 1, on the unsteady aerodynamic force reduced-order model structural nature
Norm type.
In an embodiment of the invention, when the unsteady aerodynamic force reduced-order model structure is k-th in 10 formula
The linear combination of the preceding na output and nb input nearest from kth time step is regarded in the response of spacer step as.
In an embodiment of the invention, in step 210, the unsteady aerodynamic force reduced-order model is converted into state
The process of space equation is as follows：
Make modal displacement ξ=u, mode aerodynamic coefficient fa=y, is defined as follows state vector：
xa(k)=[fa(k-1),…,fa(k-na),ξ(k-1),…,ξ(k-nb+1)]T 11
Then 10 formula of formula can be written as：
In formula：
Then 12 formulas are converted into the state space equation of the continuous system of 13 formula formulas：
In an embodiment of the invention, the generalized form of rigid body dynamic equation is：
Wherein, M is general mass matrix, and Q is broad sense aerodynamic force matrix, and ξ is generalized Modal transposed matrix；
Condition, definition status variable are coupled as with two-freedomThen dynamics of rigid bodies shown in 14 formulas
Equation is variable to be turned to：
In formula, Q represents free incoming dynamic pressure；
Using discrete space to continuous space method for transformation, the continuous system state that 15 formulas can be directly converted into 16 formulas is empty
Between form：
In an embodiment of the invention, the state space by unsteady aerodynamic force reduced-order model under continuous system
13 formula of equation can obtain currently with rigid body dynamic equation after 16 formula of state equation equation of continuous space carries out feedback link
The coupled dynamic stability analysis equation of aircraft, equation are as follows：
The present invention proposes a kind of will be based on efficiently for the full flight domain multi-channel coupling dynamic stability problem analysis of aircraft
The unsteady aerodynamic force reduced-order model of adaptively sampled method solves feature square with Dynamical Equations of Rigid Body in State-Space Coupling
The method of battle array characteristic root carries out the full flight domain dynamic stability feature prediction of aircraft, can be while precision of prediction is ensured significantly
It reduces and calculates cost, so as to be easy to carry out the qualitative and quantitative study of the full flight domain dynamic stability characteristic of aircraft, obtain to flight
Device is designed with the design direction of directive significance.
Description of the drawings
Fig. 1 is the step flow chart of one embodiment of the present invention；
Fig. 2 is the analysis process schematic diagram of one embodiment of the present invention；
Fig. 3 is the displacement signal schematic diagram of one embodiment of the present invention；
Fig. 4 is the example Waverider of one embodiment of the present invention；
Fig. 5 is the adaptive the probability distribution of samples points of one embodiment of the present invention；
Fig. 6 is the horizontal course coupling dynamic stability modal distribution in full flight domain of one embodiment of the present invention.
Specific embodiment
As shown in Figure 1, 2, full flight domain aerocraft dynamic stability analysis method disclosed by the invention is complete firstly the need of establishing
Then flight domain unsteady aerodynamic force agency-reduced-order model carries out the analysis of State-Space Coupling dynamic stability for full flight domain,
Then dynamic stability feature of the aircraft in entire flight domain can be obtained.Its process includes the following steps in general manner：
Step 100, for the Aerodynamic Characteristics of required research, determine in aimed at precision requirement and calculation amount
Limit, and choose the initial sample point of predetermined quantity in full flight domain scope；
Goal required precision refers to the precision finally to be realized after interpolation processing of the model, and in calculation amount
Limit is then the number for adding candidate's sampled point.
Here the border vertices of parameter space (full flight domain scope) should be generally included in initial sample point, in order to avoid in generation
Extrapolation occurs when managing model prediction.But in the case of parameter dimensions are very high, all vertex are included in initial sample point, similary meeting
It meets with " dimension curse ".Such as the parameter space for one 12 dimension, total vertex number are 4096 (i.e. 212), this quantity is very
To being more than total number of samples, it is clear that be worthless in the application.In this case, it will usually super vertical using optimization Latin
The random devices such as side's sampling (Optimized-Latin Hypercube Sampling, OLHS), one group of generation substantially can be complete
The sample point of the whole entire parameter space of covering carries out sampling initialization.
Step 110, in initial sample point x(i)Place, the training signal sequence u (x given to one group(i), k) and it uses
CFD calculates the when ordinal series y (x of the unsteady aerodynamic force for obtaining aircraft forced movement and torque coefficient(i), k), then with u
(x(i), k) and y (x(i), k) and it is sample data, structure Kriging agent models are as initial target proxy model；
So-called training signal refers to that aircraft meets the rigid motion time-varying process of certain frequency characteristic, usually
" 3211 " signal can meet frequency characteristic requirement.For example, to predict rolling/yaw coupled dynamic stability characteristic of aircraft,
The flight force and moment of rolling/yaw coupled motions process is then carried out at the same time by CFD solution aircraft.Input signal is using figure
" 3211 " signal shown in 2, " x1 " represents rolling angular displacement in figure, and " x2 " represents yaw displacement, calculate the aerodynamic force of acquisition
It is to export with torque.
Although here by the use of aircraft as object but it is also possible to be aerofoil profile as object.
Step 120, structure is assessed initial sample point, with reference to agent model to determine present sample precision；
In order to reduce the calculation amount of unsteady CFD, except structure target proxy model (Desired Model, DM) outside,
Also one is constructed by other interpolation methods simultaneously to adopt for auxiliary with reference to agent model (Reference Model, RM)
Sample.In each round is adaptively sampled, it is assumed that RM is predicted as true value at sampling candidate point, then the prediction of RM can be used to comment
Estimate prediction errors of the DM in same point.In entire sampling process, unsteady CFD calculating is served only for training new sample point,
Without calculating test sample point, therefore it can greatly save calculating and expend, be built so as to improve target proxy model
Imitate rate.
Step 130, new candidate's sampled point is generated using test design method；And pass through target proxy model and with reference to generation
Reason model is respectively calculated and assesses to all candidate's sampled points, and target proxy model is added in obtain each candidate's sampled point
With the precision changing value generated after reference agent model；
Here adopted using existing (Design of Experiment) DoE methods come the random candidate for generating predetermined quantity
Sampling point, the quantity of specific candidate's sampled point can arbitrarily be set as needed.
In the step, each candidate's sampled point is substituted into DM respectively and is calculated, is adopted with obtaining corresponding candidate in advance
Sampling point adds in the influence value brought after model to model accuracy, and then the influence value of candidate's sampled point is assessed by RM,
It is final available using influence value as the sequence of all candidate's sampled points to put in order.
When carrying out agent model verification, there are many error metrics indexs, such as Nash-Sutcliffe efficiency metrics
(Nash-Sutcliffe efficiency measure, NSEM), is defined as：
In formula,Represent predicted value, f (xi) represent true value,Represent f (xi) average.
Step 140, the candidate for meeting present sample required precision in precision changing value is chosen by adaptively sampled criterion
Sample point of the sampled point as next addition target proxy model and with reference to agent model is trained；
Adaptively sampled criterion (Sampling criteria) in this step predicts the probability distribution of samples points and agent model
Precision important.In sampling optimization, more unified viewpoint is so that sample point can be realized " full of sky at present
Between ", and the demand of " local strengthening " can be taken into account.It is the measurement to sample point spatial distribution full of space, sample point is full of parameter
Space is conducive to improve the global prediction precision of agent model；And local strengthening is the measurement to function response, sample point is at certain
Local strengthening in a little important sampling domains, can enhance the ability that agent model predicts the nonlinear characteristic that function responds.
The adaptively sampled criterion has sampled point distance in specified model when choosing candidate's sampled point, it is necessary to reject relatively
Candidate's sampled point in enclosing and outside specified range.For the adaptive selection of sample point, it is preferred that emphasis is adjust sample point point
The geometrical constraint of cloth should avoid the new existing sampled point of candidate's sampled point distance too near, prevent that its is too far again, sample point exists
When some or certain several local sampling domains are gathered, agent model over-fitting can be caused, and too far when may weaken and act on behalf of mould
Type captures the ability of function nonlinearity characteristic.
In addition, when candidate's sampled point occurs with having the phenomenon that sampled point overlaps, candidate's sampled point is rejected.
Therefore, in this step, the candidate's sampled point selected by adaptively sampled criterion is not necessarily to "current" model
Precision influence maximum, but be distributed most rational.
Step 150, step 130 is repeated to 140, then whether judges current goal agent model precision in each repeat
Whether up to standard and calculation amount reaches the upper limit, and previous cycle is exited when meeting one of condition, completes current goal agency
The sampling of model and the training sample that full flight domain unsteady aerodynamic force agency-reduced-order model is determined with this；
When increasing sampled point, abovementioned steps are repeated, it can be on the basis of new sampled point be added every time, according to adaptive
Sampling criterion selects the optimum sampling point in next group candidate's sampled point.
When the precision of increased candidate's sampled point meets the required precision of setting, the supplement of sampled point can be stopped.Or
It is in the case of precision backlog demand, but increased sampled point quantity is met when requiring, and can also stop sampled point
Supplement, above-mentioned setting can be to avoid using the excessive calculating time.
After establishing training sample, for any new design point, you can obtain the input of the design point by Kriging interpolation
Export signal.
In abovementioned steps, the function measurement of specific adaptively sampled criterion is：
C(ξi)=δ (ξi)*ρ(ξi) (1)
Wherein, δ (ξi) RM and DM is characterized in ξiThe deviation of the prediction at place is known as deflection function (Deviation
Function, DF).
Method when candidate samples point updates is as follows：
1) value of adaptively sampled criterion C is assessed at the candidate's sampled point being generated in advance at all Ncand；
2) the value descending arrangement for the C that above-mentioned candidate's sampled point is obtained according to previous step, first candidate sampling therein
Point is the state point for observing maximum C values, enters to elect new sample point as, waits master mould training.
3) during each round is adaptively sampled, repeat the above steps, until there are Nuppc (Number of updated
Points per circle) until a new sample point generates.
δ(ξi) calculation it is as follows：
Wherein x(j)It is existing sample point,WithIt is the prediction of RM and DM respectively；As candidate samples
Point ξiHas sample point x with certain(j)During coincidence, δ is forced zero setting, to eliminate the influence of rounding error；
ρ(ξi) expression formula it is as follows：
Wherein, molecule is candidate's sampled point ξiThe minimum standardization Euclidean distance between all existing sample points；Denominator
It is candidate's sampled point ξiThe closest distance of expectation between all existing sample points.
ρ(ξi) as the ratio between minimum range and desired distance, it is to candidate's sampled point and existing sample point spatial distribution
The measurement of pattern can effectively distinguish non-sampling area and sampling area.After parameter area determines, desired distance rexpOnly
It is related with sample point sum Ns.Desired distance r during adaptively sampledexpIt can be continuous with the increase of sample point sum Ns
Update.In addition the ratio between minimum range and desired distance ρ (ξi) value can be adjusted by index q, the value for increasing q is enhanced and not sampled
Difference between region and sampling area.For the adaptively sampled choosing then of sample point, the meaning of this parameter is to adjust
The geometrical constraint of the probability distribution of samples points is saved, big q means the distance limitation enhancing between sample point, so as to effective
Sample point is avoided to gather in some or certain several local sampling domains, causes agent model over-fitting.However excessive q is unfavorable for sample
This local enhancement sampling, may weaken the ability that agent model captures function nonlinearity characteristic.It is it can be seen that minimum
The ratio between distance and desired distance ρ (ξi) separation of sample point spatial distribution/gather characteristic is known, therefore be referred to as separating letter
Number.
Further, the closest distance of expectation in abovementioned steps can be given by：
Wherein Ns is all existing sample point numbers in parameter space, and A is the geometry of studied parameter area standardization
Characterization, calculates as follows：
x(upper)And x(lower)The upper and lower border of each variable dimension in expression parameter space respectively：
S is Nvar×NvarDiagonal matrix is tieed up, wherein l-th of diagonal entry corresponds to l-th of variable dimension xlIt is upper following
The standard deviation on boundary：
μ in formulalIt is xlThe average of up-and-down boundary：
Step 200, a large amount of random distributions, the design full of full flight domain are generated in full flight domain with arbitrary sampling method
Point；
Step 210, to each design point, the input for obtaining each design point using Kriging interpolation in training sample is defeated
Go out signal, the discrete space that each input/output signal is formed is recycled to determine unsteady aerodynamic force reduced-order model；
Wherein y (k) represents the pneumatic force vector of broad sense of kth time step, and u (k-i) represents the broad sense position of-i time steps of kth
Move input vector；AiAnd BiFor coefficient matrix to be identified；Na and nb is respectively the delay exponent number for exporting and inputting, and characterization is pneumatic
The unsteady aerodynamic effect of power.
It samples linear least squares method and solves undetermined coefficient matrix AiAnd BiJust obtain required unsteady aerodynamic force mould
Type, rule of thumb, na and nb, which generally take, 3~5 can meet required precision.
Step 220, unsteady aerodynamic force reduced-order model is converted into the state space equation of continuous space；By rigid body dynamic
The equations turned state equation under continuous space；By the state space equation and rigid body dynamic of unsteady aerodynamic force reduced-order model
The state equation of equation carries out feedback link to get to the coupled dynamic stability analysis equation of current flight device；It is unsteady pneumatic
The conversion process of power reduced-order model is as follows：
Make modal displacement ξ=u, mode aerodynamic coefficient fa=y, is defined as follows state vector：
xa(k)=[fa(k-1),…,fa(k-na),ξ(k-1),…,ξ(k-nb+1)]T (11)
Then formula (10) can be written as：
In formula：
fa0For static aerodynamic coefficient, since static aerodynamic force does not influence the dynamic stability of aircraft, divide in dynamic stability
It can ignore in analysis；
Using discrete space to continuous space method for transformation, (12) can be directly converted into the continuous system state of (13) formula
Space form：
The equations turned conversion process of rigid body dynamic is as follows：
The generalized form of rigid body dynamic equation is：
Wherein, M is general mass matrix, and Q is broad sense aerodynamic force matrix, and ξ is generalized Modal transposed matrix；
(such as rolling/yaw coupling), definition status variable by taking two-freedom couples as an exampleThen
(14) rigid body dynamic model shown in formula is variable turns to：
In formula, Q represents free incoming dynamic pressure.
Using discrete space to continuous space method for transformation, (15) can be directly converted into the continuous system state of (16) formula
Space form：
The coupled dynamic stability analysis equation of current flight device is as follows：
Step 230, the eigenmatrix characteristic root of coupled dynamic stability analysis equation, the characterization system resistance of characteristic root real part are solved
Buddhist nun, imaginary part characterization system frequency；Wherein, when all characteristic root real parts are all negative, the aircraft dynamic stability of the design point；When
When there is the characteristic root of positive real part, the aircraft of the design point moves unstable；
Step 240, after the dynamic stability feature that the aircraft of each design point is obtained by above-mentioned steps, you can flown
Dynamic stability feature of the row device in entire flight domain.
Exemplary embodiment
Using above-mentioned aircraft coupled dynamic stability characteristic analysis method, study Waverider aircraft as shown in Figure 4 and exist
Horizontal course coupling dynamic stability characteristic in full flight domain.Flight Design parameter elects Mach number M, height H and angle of attack, Mach number as
Excursion is 5~20, and height change scope is 25~70km, and angle of attack excursion is -10 °~20 °.
First, the initialization sample points of agency-reduced-order model are N0=30, adaptively sampled other parameter is q=
2、Ns1=20, Nm=10, p=1 and Ncand=6000.Model output dimension N in this examplerep=3, to each response in every wheel
It samples and 3 sample points is updated in iteration, therefore it is N often to take turns newer sample point sumuppc=9.In order to save unsteady CFD
Calculate take, when sample point sum more than 200 or it is adaptively sampled in two respond NSEM (Nash-Sutcliffe efficiency degree
Amount) assess termination sampling when meeting or exceeding 0.9.The final training sample points of this example are Ns=30+9 × 16=174,
The middle newer sample points of adaptive approach are 144, and the probability distribution of samples points is as shown in Figure 5.
Then, using " 3211 " displacement signal shown in Fig. 3 as input signal, CFD is substituted into each training sample point and is solved
In acquire the output of corresponding with input signal aerodynamic force, including rolling moment coefficient and yawing moment coefficient.
Then, 11880 design sample points are chosen in the design space of flight domain by Monte Carlo sampling, passes through agency
Model obtains the input/output signal of each design point, is substituted into formula (10), coefficient matrix is acquired by linear least squares method
AiAnd Bi, you can determine unsteady aerodynamic force reduced-order model.
Finally, the implementation of step 210~260 is passed through, you can acquire the eigenmatrix characteristic root of each design point, also just obtain
Obtained horizontal course coupled dynamic stability feature of the aircraft in full flight domain.Horizontal course coupling dynamic stability mode is entirely being flown
Distribution in domain is as shown in fig. 6, wherein the benchmark angle of attack of (a) figure is α0=0 °, the benchmark angle of attack in (b) figure is α0=10 °, figure
In " Dutch Roll Divergence " represent Dutch roll Diverging mode, " Dutch Roll Convergence " represent Holland
Rolling convergence mode, " Static Divergence " represents rolling Diverging mode.
Above-described embodiment is to present through the full flight domain aerocraft dynamic stability analysis method of this patent proposition to divide
Analyse whole process of the aircraft in full flight domain coupling dynamic stability characteristic.
So far, although those skilled in the art will appreciate that detailed herein have shown and described the exemplary of the present invention
Embodiment still, without departing from the spirit and scope of the present invention, still can be determined directly according to the present disclosure
Or derive many other variations or modifications consistent with the principles of the invention.Therefore, the scope of the invention should be understood and defined as
Cover other all these variations or modifications.
Claims (11)
1. a kind of full flight domain aerocraft dynamic stability analysis method, which is characterized in that include the following steps：
Step 100, for the Aerodynamic Characteristics of required research, aimed at precision requirement and the calculation amount upper limit are determined, and
In full flight domain, scope chooses the initial sample point of predetermined quantity；
Step 110, in initial sample point x(i)Place, the training signal sequence u (x given to one group(i), k) and it is calculated using CFD
Obtain the unsteady aerodynamic force of aircraft forced movement and the when ordinal series y (x of torque coefficient(i), k), then with u (x(i), k) and
y(x(i), k) and it is sample data, structure Kriging agent models are as initial target proxy model；
Step 120, structure is assessed initial sample point, with reference to agent model to determine present sample precision；
Step 130, new candidate's sampled point is generated using test design method；And pass through target proxy model and refer to and act on behalf of mould
Type is respectively calculated and assesses to all candidate's sampled points, and target proxy model and ginseng are added in obtain each candidate's sampled point
Examine the precision changing value generated after agent model；
Step 140, the candidate's sampling for meeting present sample required precision in precision changing value is chosen by adaptively sampled criterion
Sample point of the point as next addition target proxy model and with reference to agent model is trained；
Step 150, step 130 is repeated to 140, then judges whether current goal agent model precision is up to standard in each repeat
Whether reach the upper limit with calculation amount, and previous cycle is exited when meeting one of condition, complete current goal agent model
Sampling and the training sample of full flight domain unsteady aerodynamic force agency-reduced-order model is determined with this；
Step 200, a large amount of random distributions, the design point full of full flight domain are generated in full flight domain with arbitrary sampling method；
Step 210, to each design point, the input and output for being obtained each design point using Kriging interpolation in training sample are believed
Number, the discrete space that each input/output signal is formed is recycled to determine unsteady aerodynamic force reduced-order model；
Step 220, unsteady aerodynamic force reduced-order model is converted into the state space equation of continuous space；By rigid body dynamic equation
The state equation being converted under continuous space；By the state space equation of unsteady aerodynamic force reduced-order model and rigid body dynamic equation
State equation carry out feedback link to get to the coupled dynamic stability analysis equation of current flight device；
Step 230, the eigenmatrix characteristic root of coupled dynamic stability analysis equation is solved, characteristic root real part characterizes system damping,
Imaginary part characterizes system frequency；Wherein, when all characteristic root real parts are all negative, the aircraft dynamic stability of the design point；When going out
Now during the characteristic root of positive real part, the aircraft of the design point moves unstable；
Step 240, after the dynamic stability feature that the aircraft of each design point is obtained by above-mentioned steps, you can obtain aircraft
Dynamic stability feature in entire flight domain.
2. full flight domain aerocraft dynamic stability analysis method according to claim 1, which is characterized in that
In the step 100, initial sample point is generated using OLHS random devices, and the complete covering of needs is entire described complete winged
Row domain space and its border vertices.
3. full flight domain aerocraft dynamic stability analysis method according to claim 1, which is characterized in that
In the step 140, the adaptively sampled criterion is when choosing candidate's sampled point, it is necessary to reject relatively existing sampling
Point distance within the specified range with candidate's sampled point outside specified range.
4. full flight domain aerocraft dynamic stability analysis method according to claim 3, which is characterized in that
The function of the adaptively sampled criterion is measured：
C(ξi)=δ (ξi)*ρ(ξi) 1
Wherein, δ (ξi) characterize with reference to agent model and target proxy model in ξiThe deviation of the prediction at place；
δ(ξi) calculation it is as follows：
<mrow>
<mi>&delta;</mi>
<mrow>
<mo>(</mo>
<msub>
<mi>&xi;</mi>
<mi>i</mi>
</msub>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mfenced open = "{" close = "">
<mtable>
<mtr>
<mtd>
<mrow>
<mo>|</mo>
<msub>
<mover>
<mi>f</mi>
<mo>^</mo>
</mover>
<mrow>
<mi>R</mi>
<mi>M</mi>
</mrow>
</msub>
<mrow>
<mo>(</mo>
<msub>
<mi>&xi;</mi>
<mi>i</mi>
</msub>
<mo>)</mo>
</mrow>
<mo>-</mo>
<msub>
<mover>
<mi>f</mi>
<mo>^</mo>
</mover>
<mrow>
<mi>D</mi>
<mi>M</mi>
</mrow>
</msub>
<mrow>
<mo>(</mo>
<msub>
<mi>&xi;</mi>
<mi>i</mi>
</msub>
<mo>)</mo>
</mrow>
<mo>|</mo>
<mo>,</mo>
</mrow>
</mtd>
<mtd>
<mrow>
<msub>
<mi>&xi;</mi>
<mi>i</mi>
</msub>
<mo>&NotEqual;</mo>
<msup>
<mi>x</mi>
<mrow>
<mo>(</mo>
<mi>j</mi>
<mo>)</mo>
</mrow>
</msup>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mn>0</mn>
<mo>,</mo>
</mrow>
</mtd>
<mtd>
<mrow>
<msub>
<mi>&xi;</mi>
<mi>i</mi>
</msub>
<mo>=</mo>
<msup>
<mi>x</mi>
<mrow>
<mo>(</mo>
<mi>j</mi>
<mo>)</mo>
</mrow>
</msup>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>,</mo>
<mi>j</mi>
<mo>=</mo>
<mn>1</mn>
<mo>,</mo>
<mn>2</mn>
<mo>,</mo>
<mn>...</mn>
<mo>,</mo>
<msub>
<mi>N</mi>
<mi>s</mi>
</msub>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mn>2</mn>
</mrow>
Wherein x(j)It is existing sample point,WithIt is with reference to agent model and target proxy model respectively
Prediction；As candidate's sampled point ξiHas sample point x with certain(j)During coincidence, δ is forced zero setting, to eliminate the shadow of rounding error
It rings；
ρ(ξi) expression formula it is as follows：
<mrow>
<mi>&rho;</mi>
<mrow>
<mo>(</mo>
<msub>
<mi>&xi;</mi>
<mi>i</mi>
</msub>
<mo>)</mo>
</mrow>
<mo>=</mo>
<msup>
<mrow>
<mo>(</mo>
<mfrac>
<mrow>
<msub>
<mi>r</mi>
<mi>min</mi>
</msub>
<mrow>
<mo>(</mo>
<msub>
<mi>&xi;</mi>
<mi>i</mi>
</msub>
<mo>,</mo>
<msup>
<mi>x</mi>
<mrow>
<mo>(</mo>
<mi>j</mi>
<mo>)</mo>
</mrow>
</msup>
<mo>)</mo>
</mrow>
</mrow>
<mrow>
<msub>
<mi>r</mi>
<mi>exp</mi>
</msub>
<mrow>
<mo>(</mo>
<msub>
<mi>&xi;</mi>
<mi>i</mi>
</msub>
<mo>,</mo>
<msup>
<mi>x</mi>
<mrow>
<mo>(</mo>
<mi>j</mi>
<mo>)</mo>
</mrow>
</msup>
<mo>)</mo>
</mrow>
</mrow>
</mfrac>
<mo>)</mo>
</mrow>
<mi>q</mi>
</msup>
<mo>,</mo>
<mi>q</mi>
<mo>></mo>
<mn>0</mn>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mn>3</mn>
</mrow>
Wherein, molecule is candidate's sampled point ξiThe minimum standardization Euclidean distance between all existing sample points；Denominator is candidate
Sampled point ξiThe closest distance of expectation between all existing sample points.
5. full flight domain aerocraft dynamic stability analysis method according to claim 4, which is characterized in that
It is described it is expected that closest distance is given by：
<mrow>
<msub>
<mi>r</mi>
<mi>exp</mi>
</msub>
<mrow>
<mo>(</mo>
<msub>
<mi>&xi;</mi>
<mi>i</mi>
</msub>
<mo>,</mo>
<msup>
<mi>x</mi>
<mrow>
<mo>(</mo>
<mi>j</mi>
<mo>)</mo>
</mrow>
</msup>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mfrac>
<mn>1</mn>
<mrow>
<mn>2</mn>
<msqrt>
<mrow>
<mi>N</mi>
<mi>s</mi>
<mo>/</mo>
<mi>A</mi>
</mrow>
</msqrt>
</mrow>
</mfrac>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mn>4</mn>
</mrow>
Wherein Ns is all existing sample point numbers in full flight domain, and A is the geometry characterization of studied parameter area standardization,
It calculates as follows：
<mrow>
<mi>A</mi>
<mo>=</mo>
<msqrt>
<mrow>
<mo>(</mo>
<msup>
<mi>x</mi>
<mrow>
<mo>(</mo>
<mi>u</mi>
<mi>p</mi>
<mi>p</mi>
<mi>e</mi>
<mi>r</mi>
<mo>)</mo>
</mrow>
</msup>
<mo>-</mo>
<msup>
<mi>x</mi>
<mrow>
<mo>(</mo>
<mi>o</mi>
<mi>w</mi>
<mi>e</mi>
<mi>r</mi>
<mo>)</mo>
</mrow>
</msup>
<mo>)</mo>
<mo>*</mo>
<msup>
<mi>S</mi>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msup>
<mo>*</mo>
<msup>
<mrow>
<mo>(</mo>
<msup>
<mi>x</mi>
<mrow>
<mo>(</mo>
<mi>u</mi>
<mi>p</mi>
<mi>p</mi>
<mi>e</mi>
<mi>r</mi>
<mo>)</mo>
</mrow>
</msup>
<mo>-</mo>
<msup>
<mi>x</mi>
<mrow>
<mo>(</mo>
<mi>l</mi>
<mi>o</mi>
<mi>w</mi>
<mi>e</mi>
<mi>r</mi>
<mo>)</mo>
</mrow>
</msup>
<mo>)</mo>
</mrow>
<mi>T</mi>
</msup>
</mrow>
</msqrt>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mn>5</mn>
</mrow>
x(upper)And x(lower)The upper and lower border of each variable dimension in expression parameter space respectively：
<mrow>
<msup>
<mi>x</mi>
<mrow>
<mo>(</mo>
<mi>u</mi>
<mi>p</mi>
<mi>p</mi>
<mi>e</mi>
<mi>r</mi>
<mo>)</mo>
</mrow>
</msup>
<mo>=</mo>
<mo>&lsqb;</mo>
<msubsup>
<mi>x</mi>
<mn>1</mn>
<mrow>
<mo>(</mo>
<mi>u</mi>
<mi>p</mi>
<mi>p</mi>
<mi>e</mi>
<mi>r</mi>
<mo>)</mo>
</mrow>
</msubsup>
<mo>,</mo>
<msubsup>
<mi>x</mi>
<mn>2</mn>
<mrow>
<mo>(</mo>
<mi>u</mi>
<mi>p</mi>
<mi>p</mi>
<mi>e</mi>
<mi>r</mi>
<mo>)</mo>
</mrow>
</msubsup>
<mo>,</mo>
<mo>...</mo>
<mo>,</mo>
<msubsup>
<mi>x</mi>
<msub>
<mi>N</mi>
<mi>var</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>u</mi>
<mi>p</mi>
<mi>p</mi>
<mi>e</mi>
<mi>r</mi>
<mo>)</mo>
</mrow>
</msubsup>
<mo>&rsqb;</mo>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mn>6</mn>
</mrow>
<mrow>
<msup>
<mi>x</mi>
<mrow>
<mo>(</mo>
<mi>l</mi>
<mi>o</mi>
<mi>w</mi>
<mi>e</mi>
<mi>r</mi>
<mo>)</mo>
</mrow>
</msup>
<mo>=</mo>
<mo>&lsqb;</mo>
<msubsup>
<mi>x</mi>
<mn>1</mn>
<mrow>
<mo>(</mo>
<mi>l</mi>
<mi>o</mi>
<mi>w</mi>
<mi>e</mi>
<mi>r</mi>
<mo>)</mo>
</mrow>
</msubsup>
<mo>,</mo>
<msubsup>
<mi>x</mi>
<mn>2</mn>
<mrow>
<mo>(</mo>
<mi>l</mi>
<mi>o</mi>
<mi>w</mi>
<mi>e</mi>
<mi>r</mi>
<mo>)</mo>
</mrow>
</msubsup>
<mo>,</mo>
<mo>...</mo>
<mo>,</mo>
<msubsup>
<mi>x</mi>
<msub>
<mi>N</mi>
<mi>var</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>l</mi>
<mi>o</mi>
<mi>w</mi>
<mi>e</mi>
<mi>r</mi>
<mo>)</mo>
</mrow>
</msubsup>
<mo>&rsqb;</mo>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mn>7</mn>
</mrow>
S is Nvar×NvarDiagonal matrix is tieed up, wherein l-th of diagonal entry corresponds to l-th of variable dimension xlThe mark of up-and-down boundary
It is accurate poor：
<mrow>
<msub>
<mi>s</mi>
<mi>l</mi>
</msub>
<mo>=</mo>
<msqrt>
<mrow>
<mfrac>
<mn>1</mn>
<mrow>
<mn>2</mn>
<mo>-</mo>
<mn>1</mn>
</mrow>
</mfrac>
<mo>&lsqb;</mo>
<msup>
<mrow>
<mo>(</mo>
<msubsup>
<mi>x</mi>
<mi>l</mi>
<mrow>
<mo>(</mo>
<mi>u</mi>
<mi>p</mi>
<mi>p</mi>
<mi>e</mi>
<mi>r</mi>
<mo>)</mo>
</mrow>
</msubsup>
<mo>-</mo>
<msub>
<mi>&mu;</mi>
<mi>l</mi>
</msub>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mo>+</mo>
<msup>
<mrow>
<mo>(</mo>
<msubsup>
<mi>x</mi>
<mi>l</mi>
<mrow>
<mo>(</mo>
<mi>l</mi>
<mi>o</mi>
<mi>w</mi>
<mi>e</mi>
<mi>r</mi>
<mo>)</mo>
</mrow>
</msubsup>
<mo>-</mo>
<msub>
<mi>&mu;</mi>
<mi>l</mi>
</msub>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mo>&rsqb;</mo>
</mrow>
</msqrt>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mn>8</mn>
</mrow>
μ in formulalIt is xlThe average of up-and-down boundary：
<mrow>
<msub>
<mi>&mu;</mi>
<mi>l</mi>
</msub>
<mo>=</mo>
<mfrac>
<mn>1</mn>
<mn>2</mn>
</mfrac>
<mrow>
<mo>(</mo>
<msubsup>
<mi>x</mi>
<mi>l</mi>
<mrow>
<mo>(</mo>
<mi>u</mi>
<mi>p</mi>
<mi>p</mi>
<mi>e</mi>
<mi>r</mi>
<mo>)</mo>
</mrow>
</msubsup>
<mo>-</mo>
<msubsup>
<mi>x</mi>
<mi>l</mi>
<mrow>
<mo>(</mo>
<mi>l</mi>
<mi>o</mi>
<mi>w</mi>
<mi>e</mi>
<mi>r</mi>
<mo>)</mo>
</mrow>
</msubsup>
<mo>)</mo>
</mrow>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mn>9</mn>
</mrow>
6. full flight domain aerocraft dynamic stability analysis method according to claim 1, which is characterized in that
Unsteady aerodynamic force reduced-order model structure in the step 210 is：
<mrow>
<mi>y</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>i</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mrow>
<mi>n</mi>
<mi>a</mi>
</mrow>
</munderover>
<msub>
<mi>A</mi>
<mi>i</mi>
</msub>
<mi>y</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>-</mo>
<mi>i</mi>
<mo>)</mo>
</mrow>
<mo>+</mo>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>i</mi>
<mo>=</mo>
<mn>0</mn>
</mrow>
<mrow>
<mi>n</mi>
<mi>b</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</munderover>
<msub>
<mi>B</mi>
<mi>i</mi>
</msub>
<mi>u</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>-</mo>
<mi>i</mi>
<mo>)</mo>
</mrow>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mn>10</mn>
</mrow>
Y (k) represents the pneumatic force vector of broad sense of kth time step in formula, and u (k-i) represents that the generalized displacement of-i time steps of kth is defeated
Enter vector；AiAnd BiFor coefficient matrix to be identified；Na and nb is respectively the delay exponent number for exporting and inputting, and characterizes aerodynamic force
Unsteady aerodynamic effect.
7. the full flight domain aerocraft dynamic stability analysis method described in as requested 6, which is characterized in that
It is a kind of steady model on the unsteady aerodynamic force reduced-order model structural nature as na=0 and nb=1；Work as na=0
It is a kind of quasi-steady model on the unsteady aerodynamic force reduced-order model structural nature and during nb ＞ 1；As na ＞ 0 and nb ＞ 1
When, it is a kind of unsteady model on the unsteady aerodynamic force reduced-order model structural nature.
8. the full flight domain aerocraft dynamic stability analysis method described in as requested 7, which is characterized in that
The unsteady aerodynamic force reduced-order model structure is that the response of k-th of time step in 10 formulas is regarded as from kth time step most
Near preceding na output and the linear combination of nb input.
9. the full flight domain aerocraft dynamic stability analysis method described in as requested 8, which is characterized in that
In step 210, the process that the unsteady aerodynamic force reduced-order model is converted into state space equation is as follows：
Make modal displacement ξ=u, mode aerodynamic coefficient fa=y, is defined as follows state vector：
xa(k)=[fa(k-1),…,fa(k-na),ξ(k-1),…,ξ(k-nb+1)]T 11
Then 10 formula of formula can be written as：
<mrow>
<mfenced open = "" close = "}">
<mtable>
<mtr>
<mtd>
<mrow>
<msub>
<mi>x</mi>
<mi>a</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
<mo>=</mo>
<msub>
<mover>
<mi>A</mi>
<mo>~</mo>
</mover>
<mi>a</mi>
</msub>
<msub>
<mi>x</mi>
<mi>a</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>+</mo>
<msub>
<mover>
<mi>B</mi>
<mo>~</mo>
</mover>
<mi>a</mi>
</msub>
<mi>&xi;</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<msub>
<mi>f</mi>
<mi>a</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<msub>
<mover>
<mi>C</mi>
<mo>~</mo>
</mover>
<mi>a</mi>
</msub>
<msub>
<mi>x</mi>
<mi>a</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>+</mo>
<msub>
<mover>
<mi>D</mi>
<mo>~</mo>
</mover>
<mi>a</mi>
</msub>
<mi>&xi;</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>+</mo>
<msub>
<mi>f</mi>
<mrow>
<mi>a</mi>
<mn>0</mn>
</mrow>
</msub>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mn>12</mn>
</mrow>
In formula：
<mrow>
<msub>
<mover>
<mi>B</mi>
<mo>~</mo>
</mover>
<mi>a</mi>
</msub>
<mo>=</mo>
<msup>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<msub>
<mi>B</mi>
<mn>0</mn>
</msub>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<mn>...</mn>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<mi>I</mi>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<mn>...</mn>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
</mtr>
</mtable>
</mfenced>
<mi>T</mi>
</msup>
<mo>,</mo>
</mrow>
<mrow>
<msub>
<mover>
<mi>C</mi>
<mo>~</mo>
</mover>
<mi>a</mi>
</msub>
<mo>=</mo>
<mo>&lsqb;</mo>
<mtable>
<mtr>
<mtd>
<msub>
<mi>A</mi>
<mn>1</mn>
</msub>
</mtd>
<mtd>
<msub>
<mi>A</mi>
<mn>2</mn>
</msub>
</mtd>
<mtd>
<mn>...</mn>
</mtd>
<mtd>
<msub>
<mi>A</mi>
<mrow>
<mi>n</mi>
<mi>a</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msub>
</mtd>
<mtd>
<msub>
<mi>A</mi>
<mrow>
<mi>n</mi>
<mi>a</mi>
</mrow>
</msub>
</mtd>
<mtd>
<msub>
<mi>B</mi>
<mn>1</mn>
</msub>
</mtd>
<mtd>
<msub>
<mi>B</mi>
<mn>2</mn>
</msub>
</mtd>
<mtd>
<mn>...</mn>
</mtd>
<mtd>
<msub>
<mi>B</mi>
<mrow>
<mi>n</mi>
<mi>b</mi>
<mo>-</mo>
<mn>2</mn>
</mrow>
</msub>
</mtd>
<mtd>
<msub>
<mi>B</mi>
<mrow>
<mi>n</mi>
<mi>b</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msub>
</mtd>
</mtr>
</mtable>
<mo>&rsqb;</mo>
<mo>,</mo>
</mrow>
<mrow>
<msub>
<mover>
<mi>D</mi>
<mo>~</mo>
</mover>
<mi>a</mi>
</msub>
<mo>=</mo>
<mo>&lsqb;</mo>
<msub>
<mi>B</mi>
<mn>0</mn>
</msub>
<mo>&rsqb;</mo>
<mo>,</mo>
</mrow>
Then 12 formulas are converted into the state space equation of the continuous system of 13 formula formulas：
<mrow>
<mfenced open = "" close = "}">
<mtable>
<mtr>
<mtd>
<mrow>
<msub>
<mover>
<mi>x</mi>
<mo>&CenterDot;</mo>
</mover>
<mi>a</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>t</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<msub>
<mi>A</mi>
<mi>a</mi>
</msub>
<msub>
<mi>x</mi>
<mi>a</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>t</mi>
<mo>)</mo>
</mrow>
<mo>+</mo>
<msub>
<mi>B</mi>
<mi>a</mi>
</msub>
<mi>&xi;</mi>
<mrow>
<mo>(</mo>
<mi>t</mi>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<msub>
<mi>f</mi>
<mi>a</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>t</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<msub>
<mi>C</mi>
<mi>a</mi>
</msub>
<msub>
<mi>x</mi>
<mi>a</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>t</mi>
<mo>)</mo>
</mrow>
<mo>+</mo>
<msub>
<mi>D</mi>
<mi>a</mi>
</msub>
<mi>&xi;</mi>
<mrow>
<mo>(</mo>
<mi>t</mi>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mn>13.</mn>
</mrow>
10. the full flight domain aerocraft dynamic stability analysis method described in as requested 9, which is characterized in that
The generalized form of rigid body dynamic equation is：
<mrow>
<mi>M</mi>
<mover>
<mi>&xi;</mi>
<mo>&CenterDot;&CenterDot;</mo>
</mover>
<mo>=</mo>
<mi>Q</mi>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mn>14</mn>
</mrow>
Wherein, M is general mass matrix, and Q is broad sense aerodynamic force matrix, and ξ is generalized Modal transposed matrix；
Condition, definition status variable are coupled as with two-freedomThen Dynamical Equations of Rigid Body shown in 14 formulas
It is variable to turn to：
<mrow>
<mfenced open = "" close = "}">
<mtable>
<mtr>
<mtd>
<mrow>
<msub>
<mover>
<mi>x</mi>
<mo>&CenterDot;</mo>
</mover>
<mi>s</mi>
</msub>
<mo>=</mo>
<msub>
<mi>A</mi>
<mi>s</mi>
</msub>
<msub>
<mi>x</mi>
<mi>s</mi>
</msub>
<mo>+</mo>
<msub>
<mi>B</mi>
<mi>s</mi>
</msub>
<msub>
<mi>u</mi>
<mi>s</mi>
</msub>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mi>&xi;</mi>
<mo>=</mo>
<msub>
<mi>C</mi>
<mi>s</mi>
</msub>
<msub>
<mi>x</mi>
<mi>s</mi>
</msub>
<mo>+</mo>
<msub>
<mi>D</mi>
<mi>s</mi>
</msub>
<msub>
<mi>u</mi>
<mi>s</mi>
</msub>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mn>15</mn>
</mrow>
In formula, Q represents free incoming dynamic pressure；
Using discrete space to continuous space method for transformation, 15 formulas can be directly converted into the continuous system state space shape of 16 formulas
Formula：
<mrow>
<mfenced open = "" close = "}">
<mtable>
<mtr>
<mtd>
<mrow>
<msub>
<mover>
<mi>x</mi>
<mo>&CenterDot;</mo>
</mover>
<mi>s</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>t</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<msub>
<mi>A</mi>
<mi>s</mi>
</msub>
<msub>
<mi>x</mi>
<mi>s</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>t</mi>
<mo>)</mo>
</mrow>
<mo>+</mo>
<msub>
<mi>qB</mi>
<mi>s</mi>
</msub>
<msub>
<mi>f</mi>
<mi>a</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>t</mi>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mi>&xi;</mi>
<mrow>
<mo>(</mo>
<mi>t</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<msub>
<mi>C</mi>
<mi>s</mi>
</msub>
<msub>
<mi>x</mi>
<mi>s</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>t</mi>
<mo>)</mo>
</mrow>
<mo>+</mo>
<msub>
<mi>qD</mi>
<mi>s</mi>
</msub>
<msub>
<mi>f</mi>
<mi>a</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>t</mi>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mn>16.</mn>
</mrow>
11. the full flight domain aerocraft dynamic stability analysis method described in as requested 10, which is characterized in that
By state space equation 13 formula of the unsteady aerodynamic force reduced-order model under continuous system, with rigid body dynamic equation continuous
16 formula of state equation equation in space is the coupled dynamic stability analysis equation that can obtain current flight device after carrying out feedback link,
Equation is as follows：
<mrow>
<mfenced open = "{" close = "}">
<mtable>
<mtr>
<mtd>
<msub>
<mover>
<mi>x</mi>
<mo>&CenterDot;</mo>
</mover>
<mi>s</mi>
</msub>
<mo>(</mo>
<mi>t</mi>
<mo>)</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<msub>
<mover>
<mi>x</mi>
<mo>&CenterDot;</mo>
</mover>
<mi>a</mi>
</msub>
<mo>(</mo>
<mi>t</mi>
<mo>)</mo>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>=</mo>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mrow>
<msub>
<mi>A</mi>
<mi>s</mi>
</msub>
<mo>+</mo>
<mi>q</mi>
<mo>&CenterDot;</mo>
<msub>
<mi>B</mi>
<mi>s</mi>
</msub>
<msub>
<mi>D</mi>
<mi>a</mi>
</msub>
<msub>
<mi>C</mi>
<mi>s</mi>
</msub>
</mrow>
</mtd>
<mtd>
<mrow>
<mi>q</mi>
<mo>&CenterDot;</mo>
<msub>
<mi>B</mi>
<mi>s</mi>
</msub>
<msub>
<mi>C</mi>
<mi>a</mi>
</msub>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<msub>
<mi>B</mi>
<mi>a</mi>
</msub>
<msub>
<mi>C</mi>
<mi>s</mi>
</msub>
</mrow>
</mtd>
<mtd>
<msub>
<mi>A</mi>
<mi>a</mi>
</msub>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>&CenterDot;</mo>
<mfenced open = "{" close = "}">
<mtable>
<mtr>
<mtd>
<msub>
<mi>x</mi>
<mi>s</mi>
</msub>
<mo>(</mo>
<mi>t</mi>
<mo>)</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<msub>
<mi>x</mi>
<mi>a</mi>
</msub>
<mo>(</mo>
<mi>t</mi>
<mo>)</mo>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mn>17.</mn>
</mrow>
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 true CN108121856A (en) | 2018-06-05 |
| CN108121856B 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) |
Cited By (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109634309A (en) * | 2019-02-21 | 2019-04-16 | 南京晓庄学院 | A kind of aircraft automatic obstacle avoiding system, method and aircraft |
| CN109753690A (en) * | 2018-12-10 | 2019-05-14 | 西北工业大学 | Nonlinear unsteady aerodynamics order reducing method based on Fluid Mechanics Computation |
| CN111159942A (en) * | 2019-12-26 | 2020-05-15 | 北京电子工程总体研究所 | Method for calculating roll damping torque of winged aircraft based on steady simulation |
| CN111338364A (en) * | 2019-11-21 | 2020-06-26 | 浙江大学 | High-precision controller for optimizing trajectory of hypersonic aerocraft with quick response |
| CN112699624A (en) * | 2021-03-24 | 2021-04-23 | 南京信息工程大学 | Trajectory calculation method under severe meteorological conditions |
| CN112711815A (en) * | 2021-03-25 | 2021-04-27 | 中国科学院自动化研究所 | Aircraft modeling and model characteristic analysis system |
| CN113609600A (en) * | 2021-10-11 | 2021-11-05 | 中国空气动力研究与发展中心计算空气动力研究所 | Multi-body separation compatibility measurement and characterization method suitable for aircraft |
| CN114323552A (en) * | 2021-11-18 | 2022-04-12 | 厦门大学 | Method for judging stability of water entering and exiting from cross-medium navigation body |
Citations (2)
| 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 |
-
2017
- 2017-12-06 CN CN201711277166.8A patent/CN108121856B/en active Active
Patent Citations (2)
| 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)
| Title |
|---|
| A. DA RONCH ET AL.: "On the generation of flight dynamics aerodynamic tables by computational fluid dynamics", 《PROGRESS IN AEROSPACE SCIENCES》 * |
| 聂雪媛 等: "基于 Kriging 代理模型的飞行器结构刚度气动优化设计", 《气体物理》 * |
| 龙腾 等: "基于计算试验设计与代理模型的飞行器近似优化策略探讨", 《机械工程学报》 * |
Cited By (12)
| 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 |
| CN109634309A (en) * | 2019-02-21 | 2019-04-16 | 南京晓庄学院 | A kind of aircraft automatic obstacle avoiding system, method and aircraft |
| CN109634309B (en) * | 2019-02-21 | 2024-03-26 | 南京晓庄学院 | Autonomous obstacle avoidance system and method for aircraft and aircraft |
| CN111338364A (en) * | 2019-11-21 | 2020-06-26 | 浙江大学 | High-precision controller for optimizing trajectory of hypersonic aerocraft with quick response |
| CN111159942A (en) * | 2019-12-26 | 2020-05-15 | 北京电子工程总体研究所 | Method for calculating roll damping torque of winged aircraft based on steady simulation |
| CN111159942B (en) * | 2019-12-26 | 2023-09-15 | 北京电子工程总体研究所 | Method for calculating rolling damping moment of winged aircraft based on steady simulation |
| CN112699624A (en) * | 2021-03-24 | 2021-04-23 | 南京信息工程大学 | Trajectory calculation method under severe meteorological conditions |
| CN112711815A (en) * | 2021-03-25 | 2021-04-27 | 中国科学院自动化研究所 | Aircraft modeling and model characteristic analysis system |
| CN112711815B (en) * | 2021-03-25 | 2021-06-25 | 中国科学院自动化研究所 | Aircraft modeling and model characteristic analysis system |
| CN113609600A (en) * | 2021-10-11 | 2021-11-05 | 中国空气动力研究与发展中心计算空气动力研究所 | Multi-body separation compatibility measurement and characterization method suitable for aircraft |
| CN113609600B (en) * | 2021-10-11 | 2021-12-14 | 中国空气动力研究与发展中心计算空气动力研究所 | Multi-body separation compatibility measurement and characterization method suitable for aircraft |
| CN114323552A (en) * | 2021-11-18 | 2022-04-12 | 厦门大学 | Method for judging stability of water entering and exiting from cross-medium navigation body |
Also Published As
| Publication number | Publication date |
|---|---|
| CN108121856B (en) | 2020-08-04 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN108121856A (en) | A kind of full flight domain aerocraft dynamic stability analysis method | |
| CN105843073B (en) | A kind of wing structure aeroelastic stability analysis method not knowing depression of order based on aerodynamic force | |
| CN107133406B (en) | Rapid search method for static voltage stability domain boundary of power system | |
| CN112016167B (en) | Aircraft aerodynamic shape design method and system based on simulation and optimization coupling | |
| CN104881510B (en) | A kind of lifting airscrew/tail-rotor aerodynamic interference numerical value emulation method | |
| CN107391891A (en) | A kind of high aspect ratio wing Optimization Design based on Model Fusion method | |
| CN106055791B (en) | Aircraft overall situation Aerodynamic optimization method based on Predictor-Correcting Algorithm | |
| CN106354695A (en) | Output-only linear time-varying structure modal parameter identification method | |
| CN114065662B (en) | Method suitable for rapidly predicting airfoil flow field with variable grid topology | |
| CN104834772B (en) | Aircraft wing based on artificial neural network/wing inverse design method | |
| AU2021103515A4 (en) | Data assimilation method based on multi-layer perceptron | |
| CN109840349B (en) | Fixed-wing aircraft gust response modeling analysis method | |
| CN105718634A (en) | Airfoil robust optimization design method based on non-probability interval analysis model | |
| CN107976908A (en) | A kind of aircraft coupled dynamic stability characteristic analysis method | |
| CN112947534A (en) | Adaptive pseudo-spectral method trajectory optimization method for depression section of hypersonic aircraft | |
| CN108038292A (en) | A kind of efficient self-adapted method of sampling based on dual-proxy technology | |
| CN112016253A (en) | High-fidelity chaotic polynomial correction method suitable for CFD uncertainty quantification | |
| CN103310060A (en) | Transonic limit cycle flutter analysis method | |
| CN111581784A (en) | Flapping wing motion parameter optimization method based on data-driven self-adaptive quasi-steady-state model | |
| CN112434380A (en) | Aircraft wing weight estimation method | |
| CN115146538A (en) | Power system state estimation method based on message passing graph neural network | |
| CN111199105A (en) | Flapping wing motion parameter optimization method | |
| CN117171894B (en) | Aircraft layout pneumatic optimization design method considering static margin constraint | |
| CN114611420A (en) | Unsteady aerodynamic force calculation precision evaluation and correction method | |
| CN103065015A (en) | Internal force path geometrical morphology based low-carbon material-saving bearing structure design method |
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 |