CN113221350B - Hypersonic aircraft transition prediction method based on global stability analysis - Google Patents

Hypersonic aircraft transition prediction method based on global stability analysis Download PDF

Info

Publication number
CN113221350B
CN113221350B CN202110507168.1A CN202110507168A CN113221350B CN 113221350 B CN113221350 B CN 113221350B CN 202110507168 A CN202110507168 A CN 202110507168A CN 113221350 B CN113221350 B CN 113221350B
Authority
CN
China
Prior art keywords
equation
flow direction
flow
curve
disturbance
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
CN202110507168.1A
Other languages
Chinese (zh)
Other versions
CN113221350A (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.)
Tianjin University
Original Assignee
Tianjin University
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 Tianjin University filed Critical Tianjin University
Priority to CN202110507168.1A priority Critical patent/CN113221350B/en
Publication of CN113221350A publication Critical patent/CN113221350A/en
Application granted granted Critical
Publication of CN113221350B publication Critical patent/CN113221350B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00Computer-aided design [CAD]
    • G06F30/20Design optimisation, verification or simulation
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00Three dimensional [3D] modelling, e.g. data description of 3D objects

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Geometry (AREA)
  • General Physics & Mathematics (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • General Engineering & Computer Science (AREA)
  • Computer Graphics (AREA)
  • Software Systems (AREA)
  • Aerodynamic Tests, Hydrodynamic Tests, Wind Tunnels, And Water Tanks (AREA)

Abstract

The invention relates to a hypersonic aircraft transition prediction method based on global stability analysis, which is realized by the following steps: calculating to obtain a laminar flow basic flow field on the surface of a target high supersonic aircraft to obtain a three-dimensional boundary layer; step two, performing two-dimensional global stability analysis on the three-dimensional boundary layer of the basic flow field obtained in the step one at different flow direction positions to obtain a plurality of growth rate cloud charts of unstable modes; step three, calculating the maximum growth factor of each unstable mode at each flow direction position obtained in the step two, representing the maximum growth factor by using a Nmax curve, and obtaining a Nmax _ all curve by using envelope lines of the Nmax curves of different modes again to represent the conservatively estimated maximum growth factor of the unstable mode; and step four, predicting the transition position of the target hypersonic aircraft according to a transition criterion Nt.

Description

Hypersonic aircraft transition prediction method based on global stability analysis
The technical field is as follows:
the invention relates to the technical field of transition prediction of aircrafts, in particular to a transition prediction method of a boundary layer of a hypersonic aircraft, and particularly relates to a transition prediction method of the hypersonic aircraft based on global stability analysis.
Background art:
under the hypersonic flight condition, the friction resistance and the heat flow of the laminar flow area and the turbulent flow area have a great difference, if the transition position can be accurately predicted, when a hypersonic flight vehicle thermal protection system is designed, the problem that the weight of the hypersonic flight vehicle is overlarge due to conservative design completely according to the turbulent flow state can be avoided, so that the flight distance and the maneuverability of the hypersonic flight vehicle are greatly improved, and the hypersonic flight vehicle thermal protection system has important significance for research and development of the hypersonic flight vehicle.
Currently, the frequently adopted transition prediction methods include: the method comprises a traditional eN method based on semi-experience of a linear stability theory, a method (PSE) for solving a parabolic stability equation, a Direct Numerical Simulation (DNS) and a transition mode, wherein the eN method is generally considered to be an effective method for predicting the transition position in engineering application.
The eN method is a transition prediction method based on a linear stability theory, and is traditionally directed at a one-dimensional profile. The transition prediction method predicts the transition by calculating the linear increase multiple of the unstable wave, and the transition criterion, namely the critical value of N, is usually given by depending on experiments or practical experience, so that the transition prediction method is a semi-theoretical semi-empirical transition prediction method.
The amplification factor in the evolution of the amplitude of a frequency disturbance along the flow direction may be e-exponential, i.e. e, of the amplification factor NNExpressed, the expression is:
Figure BDA0003058889600000011
the expression for the value of the amplification factor N is therefore:
Figure BDA0003058889600000012
wherein ξ0Is the starting position of integration, xi is any position downstream of the integration path, A0Is an initial position xi0The disturbance amplitude is set, A is the disturbance amplitude at the local xi position, -alphaiIs the spatial growth rate. Calculating N values under different frequencies, and then obtaining envelope curves of the N values under different frequencies, namely NmaxCurve line. For NmaxThe value can be used for proposing a transition criterion Nt according to experiments and practical experience, when N ismaxWhen the value reaches Nt, a transition is considered to occur.
The conventional eN method is directed to a one-dimensional elementary stream profile, and only the variations of the elementary stream profile in the normal direction are considered in the stability analysis. However, the basic flow field on the surface of the hypersonic aircraft mostly has an area which is changed drastically both along the normal direction and the span direction, and a two-dimensional basic flow section is arranged at the same flow direction position, so that the transition position of the area cannot be accurately predicted by the conventional eN method. Therefore, a transition prediction method based on two-dimensional global stability analysis is needed to be provided for the hypersonic aircraft and simultaneously considering the changes of the basic flow profile along the normal direction and the span direction.
The invention content is as follows:
the invention aims to provide a hypersonic aircraft transition prediction method based on global stability analysis, which considers the characteristic that a basic flow section changes violently along the normal direction and the span direction at the same time based on a two-dimensional global stability analysis method, makes up the defect that the conventional eN method only considers the change of the basic flow section along the normal direction, and simultaneously makes N of a plurality of unstable modes changemaxObtaining N by taking envelope again from the curvemax_allAnd (4) obtaining a transition criterion Nt for accurately predicting the transition position, and solving the problem of engineering application of the traditional eN method. The technical solution of the invention is as follows:
a transition prediction method of a hypersonic aircraft based on global stability analysis is realized by the following steps:
calculating to obtain a laminar flow basic flow field on the surface of a target high supersonic aircraft to obtain a three-dimensional boundary layer;
step two, performing two-dimensional global stability analysis on the three-dimensional boundary layer of the basic flow field obtained in the step one at different flow direction positions to obtain a plurality of growth rate cloud charts of unstable modes;
step three, calculating the maximum growth factor of each unstable mode obtained in the step two at each flow direction position, and using NmaxIs shown by a curve and is applied to N of different modesmaxObtaining N by taking envelope again from the curvemax_allA curve representing the maximum growth factor of a conservatively estimated unstable mode;
and step four, predicting the transition position of the target hypersonic aircraft according to a transition criterion Nt.
In the second step, the method for performing two-dimensional global stability analysis at different flow direction positions comprises the following steps: establishing a compressible dimensionless N-S equation in a matrix form by adopting a matrix-free two-dimensional global stability analysis method suitable for a three-dimensional boundary layer of a hypersonic aircraft; expressing the instantaneous quantity as the sum of the basic flow and the disturbance quantity, subtracting the equation of the basic flow from the equation of the instantaneous quantity, and omitting the high-order small quantity to obtain a linear disturbance equation; and then, expressing the disturbance into a flow direction wave form by using a local approximate parallel flow hypothesis, and finally obtaining a frequency-free two-dimensional global stability equation required by matrix-free two-dimensional global stability analysis.
Further, the method for performing two-dimensional global stability analysis at different flow direction positions is performed according to the following steps:
1) establishing a compressible dimensionless N-S equation in a matrix form under a Cartesian coordinate system:
Figure BDA0003058889600000021
wherein x, y, z respectively represent the flow direction, normal direction and spreading direction in a cartesian coordinate system, phi represents the instantaneous vector phi ═ p, u, v, w, T)TWhere ρ represents density, u represents flow velocity, v represents normal velocity, w represents spanwise velocity, T represents temperature, Γ0、A0、B0、C0、D0、V11、V22、V33、V12、V23、V13Coefficient matrices representing the linear part of the equation, FnRepresenting a non-linear portion;
2) expressing the instantaneous quantity as the sum of the basic flow and the disturbance quantity, i.e.
Figure BDA0003058889600000022
Wherein
Figure BDA0003058889600000023
Representing the amount of a steady elementary stream,
Figure BDA0003058889600000024
the density is expressed as a function of time,
Figure BDA0003058889600000025
the speed of the flow direction is indicated,
Figure BDA0003058889600000026
which is indicative of the normal velocity of the wire,
Figure BDA0003058889600000027
the speed of the span-wise direction is indicated,
Figure BDA0003058889600000031
represents the temperature; phi '═ of (ρ', u ', v', w ', T')TFor the disturbance amount, ρ ' represents density, u ' represents flow direction velocity, v ' represents normal velocity, w ' represents span direction velocity, and T ' represents temperature; instantaneous quantity phi and basic flow
Figure BDA0003058889600000032
The equation of the basic flow is subtracted from the equation of the instantaneous quantity, and the nonlinear high-order small quantity is omitted, so that the linear disturbance equation under a Cartesian coordinate system can be obtained, and the form is as follows:
Figure BDA0003058889600000033
wherein the variable φ 'represents the perturbation vector (ρ', u ', v', w ', T')T
Γ、A、B、C、D、Vxx、Vyy、Vzz、Vxy、Vyz、VxzCoefficient matrix representing linear disturbance equation, in which basic flow is contained
Figure BDA0003058889600000034
The information of (a);
3) the linear disturbance equation (2) obtains a linear disturbance equation under a curve coordinate system through coordinate transformation:
Figure BDA0003058889600000035
where ξ, η, ζ denote the flow direction, normal direction and spread direction respectively in a curved coordinate system, and Φ 'denotes a disturbance vector Φ' ((ρ ', u', v ', w', T))T,Γ、
Figure BDA0003058889600000036
D、Vξξ、Vηη
Figure BDA0003058889600000038
Vξη
Figure BDA0003058889600000039
Coefficient matrix representing linear disturbance equation in curve coordinate system, including basic flow
Figure BDA0003058889600000037
The expression of each coefficient matrix is
Figure BDA0003058889600000041
Figure BDA0003058889600000042
Figure BDA0003058889600000043
Vξξ=Vxx,
Figure BDA0003058889600000044
Figure BDA0003058889600000045
Vξη=ηyVxyzVxz,
Vξζ=ζyVzyzVxz,
Vηζ=2ηyζyVyy+2ηzζzVzz+(ηyζzzζy)Vyz
(4)
Wherein J represents Jacobi matrix determinant with expression of
Figure BDA0003058889600000046
4) Deriving a control equation for the matrix-free two-dimensional global stability method, i.e. a frequency-free two-dimensional global stability equation, using a local parallel flow assumption to set the disturbance in the form of a flow direction waveform
Figure BDA0003058889600000047
Wherein
Figure BDA0003058889600000048
A function of the shape representing the time t,
Figure BDA0003058889600000049
wherein
Figure BDA00030588896000000410
The density is expressed as a function of time,
Figure BDA00030588896000000411
the speed of the flow direction is indicated,
Figure BDA00030588896000000412
which is indicative of the normal velocity of the wire,
Figure BDA00030588896000000413
the speed of the span-wise direction is indicated,
Figure BDA00030588896000000414
denotes temperature, α denotes flow direction wave number, c.c. denotes complex conjugation; substituting the formula (6) into a linear disturbance equation (3) under a curve coordinate system, and sorting to obtain a control equation used by the matrix-free global stability method, namely a frequency-free two-dimensional global stability equation:
Figure BDA00030588896000000415
wherein gamma is,
Figure BDA0003058889600000051
Coefficient matrix representing control equation, in which elementary streams are included
Figure BDA0003058889600000052
And a real number form flow direction wave number alpha in a time mode, the expression of each coefficient matrix being
Figure BDA0003058889600000053
Further, in the second step, the method for obtaining the growth rate cloud charts of the plurality of unstable modes comprises: calculating a time mode characteristic value lambda of a two-dimensional global stability equation (7) by adopting an Arnoldi method, and obtaining a disturbed time mode growth rate omega through the relation between the characteristic value lambda and the frequency omegai(ii) a Then gain the growth rate-alpha of the space mode which is more in line with the actual physical situation through Gaster transformationiiWhich is the imaginary part of the flow direction wave number α in complex form in spatial mode), the form of the Gaster transform is: - αi=ωi/cgWherein c isgFor group velocity, the expression is
Figure BDA0003058889600000054
Finally canObtaining the growth rate-alpha of different flow direction positions of the basic flow field in the step one under various frequenciesiAnd a characteristic function, thereby obtaining a contour cloud picture of the growth rate of a plurality of unstable modes.
Further, in step two, for equation (7), the relationship between the eigenvalue λ and the frequency ω is
Figure BDA0003058889600000055
Wherein, step III Nmax_allThe curve is obtained as follows: integrating the growth rate of each most unstable mode along the flow direction based on the growth rate cloud pictures of a plurality of unstable modes under different frequencies to obtain the distribution of the amplitude amplification factor N value of each frequency disturbance along the flow direction coordinate, enveloping the maximum value of the N value of different frequency waves along the flow direction coordinate distribution by using a curve, and representing the maximum amplitude amplification factor which can be reached by each frequency disturbance at the mode by using an envelope curve, namely obtaining the NmaxA curve; then, taking the guard estimate, and combining N of different modesmaxThe curve is wrapped again to obtain the maximum amplitude amplification factors of a plurality of unstable modes, namely Nmax_allCurve line.
The determination method of the transition criterion Nt in the fourth step is as follows: the transition criterion Nt is obtained and corrected through flight experiments or ground wind tunnel experiment data, and the amplitude of disturbance with a certain frequency reaches e of the initial amplitude in the process of propagating towards the flow directionNtThe laminar flow develops into turbulent flow; n obtained in step threemax_allFinding the flow direction position corresponding to Nt on the curve, namely considering the transition starting position Xt of the hypersonic aircraft surface.
The invention provides a hypersonic aircraft transition prediction method based on global stability analysis, which has the advantages compared with the prior art that:
(1) the hypersonic aircraft surface has a region with a basic flow field which changes violently along the normal direction and the span direction, and the transition position of the region can be predicted by the method.
(2) The invention adopts a two-dimensional global stability analysis method, considers the changes of the basic flow profile along the normal direction and the span direction, has physical basis and has more reliable transition prediction result.
Description of the drawings:
FIG. 1 is a block diagram of the steps of the present invention
FIG. 2 is a schematic diagram of a model of a target hypersonic aircraft
FIG. 3 is a cloud chart of elementary stream flow direction velocity of a target hypersonic aircraft model
FIG. 4 is a cloud chart of growth rates of 6 most unstable modes in boundary layers of basic flow fields of an aircraft model
FIG. 5 shows N in the 6 most unstable modesmaxCurve
FIG. 6 shows N in the 6 most unstable modesmax_allCurve
The specific implementation mode is as follows:
the present invention will be described in detail with reference to the following embodiments and accompanying drawings. In the following description, for purposes of explanation and not limitation, specific details are set forth in order to provide a thorough understanding of the present invention. However, it will be apparent to one skilled in the art that the present invention may be practiced in other embodiments that depart from these specific details.
It should be noted that, in order to avoid obscuring the present invention with unnecessary details, only the device structures and/or processing steps that are closely related to the scheme according to the present invention are shown in the drawings, and other details that are not so relevant to the present invention are omitted.
The embodiment of the invention provides a hypersonic aircraft transition prediction method based on global stability analysis, which comprises the following steps of:
step one, calculating to obtain a laminar flow basic flow field of a target high supersonic aircraft;
in this step, the target hypersonic aircraft of this embodiment is an elliptical cone model with a ratio of major axis to minor axis of 3:1, the radius of the blunt tip is 0.95mm, the radius of the minor axis is 41mm, the flight mach number Ma is 6, the angle of attack is 0 degree, the sideslip angle is 0 degree, a basic flow field is calculated by Direct Numerical Simulation (DNS), the characteristic length used without quantization is the radius of the blunt tip, the characteristic speed is the free incoming flow speed, and the aircraft model is schematically shown in fig. 2;
in the step, the number of normal grids in the boundary layer is required to be more than 100 so as to meet the requirement of subsequent global stability analysis; the calculated flow direction velocity cloud chart of the elementary flow field is shown in fig. 3, and only the result of the 1/4 model is calculated due to symmetry;
step two, performing two-dimensional global stability analysis on the three-dimensional boundary layer of the basic flow field obtained in the step one at different flow direction positions to obtain a plurality of growth rate cloud charts of unstable modes;
in this step, the global stability analysis step is as follows:
1) starting from a compressible dimensionless N-S equation in matrix form in a cartesian coordinate system, the equation is as follows:
Figure BDA0003058889600000061
wherein x, y, z respectively represent the flow direction, normal direction and spreading direction in a cartesian coordinate system, phi represents the instantaneous vector phi ═ p, u, v, w, T)TWhere ρ represents density, u represents flow velocity, v represents normal velocity, w represents spanwise velocity, T represents temperature, Γ0、A0、B0、C0、D0、V11、V22、V33、V12、V23、V13Coefficient matrices representing the linear part of the equation, FnRepresenting the nonlinear part, the matrix elements are specifically referred to: transition mechanism and prediction of peristatic, Su rainbow, Zhang Yongming, supersonic/hypersonic boundary layer [ M]Beijing, science publishers;
2) expressing the instantaneous quantity as the sum of the basic flow and the disturbance quantity, i.e.
Figure BDA0003058889600000071
Wherein
Figure BDA0003058889600000072
Representing the amount of a steady elementary stream,
Figure BDA0003058889600000073
the density is expressed as a function of time,
Figure BDA0003058889600000074
the speed of the flow direction is indicated,
Figure BDA0003058889600000075
which is indicative of the normal velocity of the wire,
Figure BDA0003058889600000076
the speed of the span-wise direction is indicated,
Figure BDA0003058889600000077
represents the temperature; phi '═ of (ρ', u ', v', w ', T')TFor the disturbance amount, ρ ' represents density, u ' represents flow direction velocity, v ' represents normal velocity, w ' represents span direction velocity, and T ' represents temperature. Instantaneous quantity phi and basic flow
Figure BDA0003058889600000078
The equation of the basic flow is subtracted from the equation of the instantaneous quantity, and the nonlinear high-order small quantity is omitted, so that the linear disturbance equation under a Cartesian coordinate system can be obtained, and the form is as follows:
Figure BDA0003058889600000079
wherein the variable phi' represents a disturbance vector
Figure BDA00030588896000000710
Γ、A、B、C、D、Vxx、Vyy、Vzz、Vxy、Vyz、VxzCoefficient matrix representing linear disturbance equation, in which basic flow is contained
Figure BDA00030588896000000711
The coefficient matrix elements refer to: stability characteristics and disturbance evolution of Zhang Shaolong-hypersonic 2:1 elliptic cone boundary layer [ D]Tianjin university, 2016;
3) equation (2) obtains a linear disturbance equation under a curve coordinate system through coordinate transformation
Figure BDA00030588896000000712
Where ξ, η, ζ denote the flow direction, normal direction and spread direction respectively in a curved coordinate system, and Φ 'denotes a disturbance vector Φ' ((ρ ', u', v ', w', T))T,Γ、
Figure BDA00030588896000000713
D、Vξξ、Vηη
Figure BDA00030588896000000715
Vξη
Figure BDA00030588896000000716
Coefficient matrix representing linear disturbance equation in curve coordinate system, including basic flow
Figure BDA00030588896000000714
The expression of each coefficient matrix is
Figure BDA0003058889600000081
Figure BDA0003058889600000082
Figure BDA0003058889600000083
Vξξ=Vxx,
Figure BDA0003058889600000084
Figure BDA0003058889600000085
Vξη=ηyVxyzVxz,
Vξζ=ζyVzyzVxz,
Vηζ=2ηyζyVyy+2ηzζzVzz+(ηyζzzζy)Vyz
(4)
Wherein J represents Jacobi matrix determinant with expression of
Figure BDA0003058889600000086
4) A control equation used by a matrix-free two-dimensional global stability method, namely a frequency-free two-dimensional global stability equation, is obtained through derivation, and the time mode problem is solved. Using the assumption of local parallel flow, the perturbation is set to flow direction wave form
Figure BDA0003058889600000087
Wherein
Figure BDA0003058889600000088
A function of the shape representing the time t,
Figure BDA0003058889600000089
wherein
Figure BDA00030588896000000810
The density is expressed as a function of time,
Figure BDA00030588896000000811
the speed of the flow direction is indicated,
Figure BDA00030588896000000812
which is indicative of the normal velocity of the wire,
Figure BDA00030588896000000813
the speed of the span-wise direction is indicated,
Figure BDA00030588896000000814
denotes temperature, α denotes flow direction wave number, and c.c. denotes complex conjugation. The concepts and methods referred to herein are well known in the art, and reference is made specifically to: zhongheng, Zhao Tuff, flow stability [ M]National defense industry publishers, 2004.
Substituting the formula (6) into a linear disturbance equation (3) under a curve coordinate system, and sorting to obtain a control equation used by the matrix-free global stability method, namely a frequency-free two-dimensional global stability equation
Figure BDA00030588896000000815
Wherein gamma is,
Figure BDA0003058889600000091
Coefficient matrix representing control equation, in which elementary streams are included
Figure BDA0003058889600000092
And a real number form flow direction wave number alpha in a time mode, the expression of each coefficient matrix being
Figure BDA0003058889600000093
5) Calculating a time mode characteristic value lambda of a two-dimensional global stability equation (7) by adopting an Arnoldi method, and obtaining a disturbed time mode growth rate omega through the relation between the characteristic value lambda and the frequency omegai(ii) a Then gain the growth rate-alpha of the space mode which is more in line with the actual physical situation through Gaster transformationiiAn imaginary part flowing toward the wave number α in a complex form in a spatial mode); finally, the growth rate-alpha of different flow direction positions of the basic flow field in the step one at each frequency can be obtainediAnd a characteristic function, thereby obtaining a contour cloud picture of growth rates of a plurality of unstable modes;
in this step, the Arnoldi method is a method known in the art, and can be implemented by calling an open source ARPACK package, specifically referring to: R.B.Lehoucq, D.C.Sorensen, C.Yang.ARPACK Users Guide Solution of Large Scale Eigenvalue schemes by Implicitly corrected Arnoldi Methods [ J ]. Journal of Chemical Physics,2015,157(S14): 044119.
In this step, the relationship between the characteristic value λ and the frequency ω is
Figure BDA0003058889600000094
In particular, equation (7) relates to
Figure BDA0003058889600000095
Defines a linear evolution operator L (t, α) for a certain wavenumber α, which can be expressed as
Figure BDA0003058889600000096
Taking a fixed time interval T0Then there is
Figure BDA0003058889600000097
Wherein L (T)0α) is still a linear operator, then it has a characteristic modality, i.e.
Figure BDA0003058889600000098
Wherein the value of the lambda is a characteristic value,
Figure BDA0003058889600000099
to correspond to λA characteristic function. On the other hand, in equation (7), corresponding to a certain characteristic mode
Figure BDA00030588896000000910
From T-0 to T-T0Can also be described in complex form as frequency ω
Figure BDA00030588896000000911
Where the real part of ωrRepresenting frequency, imaginary part omegaiIndicating the rate of growth. For the characteristic modality, it can be described as
Figure BDA0003058889600000101
In summary, the relationship between the eigenvalue λ and the frequency ω is
Figure BDA0003058889600000102
In this step, the Gaster transform is in the form:
i=ωi/cg。 (14)
wherein c isgFor group velocity, the expression is
Figure BDA0003058889600000103
This method is well known in the art, with specific reference to: gate M.A note on the relationship between the temporal interaction and the spatial interaction distribution in the dynamic status [ J].Journal of Fluid Mechanics,1962,14(2):222-224。
In the step, for a target hypersonic aircraft model, cloud charts of growth rates of 6 most unstable modes distributed along the flow direction are obtained, and the cloud charts are respectively Y-Mode1, Y-Mode2, Y-Mode3, Z-Mode1, Z-Mode2 and Z-Mode3, and are shown in FIG. 4;
step three, calculating the maximum growth factor of each unstable mode at each flow direction position obtained in the step twoBy NmaxThe curve is shown and the envelope is again taken to obtain Nmax_allA curve representing the maximum growth factor of a conservatively estimated unstable mode;
in this step, it is considered that the perturbations of different frequencies in the boundary layer can continuously increase to the transition position, and the amplification factor of the amplitude of the small perturbation of a given frequency from the initial position to a certain position is given by equation (14):
Figure BDA0003058889600000104
in which ξ0Is the starting position of integration, xi is any position downstream of the integration path, A0Is an initial position xi0The disturbance amplitude is set, A is the disturbance amplitude at the local xi position, -alphaiIs the spatial growth rate.
In the step, based on a plurality of growth rate cloud pictures of unstable modes, the growth rate of the most unstable mode is integrated along the flow direction under different frequencies, the distribution of amplitude amplification factor N values of various frequency disturbances along the flow direction coordinate can be obtained, the maximum values of the N values of different frequency waves along the flow direction coordinate distribution are enveloped by a curve, the envelope represents the maximum amplitude amplification factor which can be reached by various frequency disturbances, and the N is obtainedmaxA curve; n of the 6 most unstable modes of the present embodimentmaxThe curves are shown in FIG. 5; to take a conservative estimate, the N ismaxThe curve is wrapped again to obtain the maximum amplitude amplification factors of a plurality of unstable modes, namely Nmax_alCurve l, N of this examplemax_allThe curves are shown in FIG. 6;
predicting the transition position of the surface of the target hypersonic aircraft according to a transition criterion Nt;
in this step, it is considered that the amplitude reaches e of the initial amplitude as long as the disturbance of a certain frequency is in the process of propagating to the flow directionNtThe laminar flow develops into turbulent flow; n obtained in step threemax_allFinding the flow direction position corresponding to Nt on the curve, namely considering as the transition starting position of the hypersonic aircraft surfaceXt;
In this step, the transition criterion Nt of this embodiment is 8.5, and is obtained by calibrating a wind tunnel experiment result and a numerical calculation result of a known hypersonic aircraft model.
In this step, N of the target hypersonic aircraft model is combined in this embodimentmax_allThe transition position of the 3:1 elliptical cone short axis region is predicted to be about X to 247 by the curve (fig. 6) and the transition criterion Nt to 8.5.
Features that are described and/or illustrated above with respect to one embodiment may be used in the same way or in a similar way in one or more other embodiments and/or in combination with or instead of the features of the other embodiments.
It should be emphasized that the term "comprises/comprising" when used herein, is taken to specify the presence of stated features, integers, steps or components but does not preclude the presence or addition of one or more other features, integers, steps, components or groups thereof.
The above devices and methods of the present invention can be implemented by hardware, or can be implemented by hardware and software. The present invention relates to a computer-readable program which, when executed by a logic section, enables the logic section to realize the above-described apparatus or constituent section, or to realize the above-described various methods or steps. The present invention also relates to a storage medium such as a hard disk, a magnetic disk, an optical disk, a DVD, a flash memory, or the like, for storing the above program.
The many features and advantages of these embodiments are apparent from the detailed specification, and thus, it is intended by the appended claims to cover all such features and advantages of these embodiments which fall within the true spirit and scope thereof. Further, since numerous modifications and changes will readily occur to those skilled in the art, it is not desired to limit the embodiments of the invention to the exact construction and operation illustrated and described, and accordingly, all suitable modifications and equivalents may be resorted to, falling within the scope thereof.
The invention has not been described in detail and is in part known to those of skill in the art.

Claims (3)

1. A transition prediction method of a hypersonic aircraft based on global stability analysis is realized by the following steps:
calculating to obtain a laminar flow basic flow field on the surface of a target high supersonic aircraft to obtain a three-dimensional boundary layer;
step two, performing two-dimensional global stability analysis on the three-dimensional boundary layer of the basic flow field obtained in the step one at different flow direction positions to obtain a plurality of growth rate cloud charts of unstable modes, wherein the method for performing the two-dimensional global stability analysis at different flow direction positions comprises the following steps: establishing a compressible dimensionless N-S equation in a matrix form by adopting a matrix-free two-dimensional global stability analysis method suitable for a three-dimensional boundary layer of a hypersonic aircraft; expressing the instantaneous quantity as the sum of the basic flow and the disturbance quantity, subtracting the equation of the basic flow from the equation of the instantaneous quantity, and omitting the high-order small quantity to obtain a linear disturbance equation; and then, expressing the disturbance into a flow direction wave form by using a local approximate parallel flow hypothesis, and finally obtaining a frequency-free two-dimensional global stability equation required by matrix-free two-dimensional global stability analysis, wherein the method of the cloud graph of the growth rate of a plurality of unstable modes comprises the following steps: calculating a time mode characteristic value lambda of a two-dimensional global stability equation by adopting an Arnoldi method, and obtaining a disturbed time mode growth rate omega through the relation between the characteristic value lambda and the frequency omegai(ii) a Then gain the growth rate-alpha of the space mode which is more in line with the actual physical situation through Gaster transformationi,αiIn the form of complex flow in the spatial mode, towards the imaginary part of the wave number α, the Gaster transform is of the form: - αi=ωi/cgWherein c isgFor group velocity, the expression is
Figure FDA0003441681810000011
Finally, the growth rate-alpha of different flow direction positions of the basic flow field in the step one at each frequency can be obtainediAnd a characteristic function, thereby obtaining a contour cloud picture of growth rates of a plurality of unstable modes;
step three, calculating the maximum growth factor of each unstable mode at each flow direction position obtained in the step twoBy NmaxIs shown by a curve and is applied to N of different modesmaxObtaining N by taking envelope again from the curvemax_allCurve showing the maximum growth factor of a conservatively estimated unstable mode, where Nmax_allThe curve is obtained as follows: integrating the growth rate of each most unstable mode along the flow direction based on the growth rate cloud pictures of a plurality of unstable modes under different frequencies to obtain the distribution of the amplitude amplification factor N value of each frequency disturbance along the flow direction coordinate, enveloping the maximum value of the N value of different frequency waves along the flow direction coordinate distribution by using a curve, and representing the maximum amplitude amplification factor which can be reached by each frequency disturbance at the mode by using an envelope curve, namely obtaining the NmaxA curve; then, taking the guard estimate, and combining N of different modesmaxThe curve is wrapped again to obtain the maximum amplitude amplification factors of a plurality of unstable modes, namely Nmax_allA curve;
fourthly, predicting the transition position of the target hypersonic aircraft according to a transition criterion Nt;
the method for analyzing the two-dimensional global stability at different flow direction positions is implemented according to the following steps:
1) establishing a compressible dimensionless N-S equation in a matrix form under a Cartesian coordinate system:
Figure FDA0003441681810000012
wherein x, y, z respectively represent the flow direction, normal direction and spreading direction in a cartesian coordinate system, phi represents the instantaneous vector phi ═ p, u, v, w, T)TWhere ρ represents density, u represents flow velocity, v represents normal velocity, w represents spanwise velocity, T represents temperature, Γ0、A0、B0、C0、D0、V11、V22、V33、V12、V23、V13Coefficient matrices representing the linear part of the equation, FnRepresenting a non-linear portion;
2) representing instantaneous quantities as base and disturbance quantitiesTo sum, i.e.
Figure FDA0003441681810000021
Wherein
Figure FDA0003441681810000022
Representing the amount of a steady elementary stream,
Figure FDA0003441681810000023
the density is expressed as a function of time,
Figure FDA0003441681810000024
the speed of the flow direction is indicated,
Figure FDA0003441681810000025
which is indicative of the normal velocity of the wire,
Figure FDA0003441681810000026
the speed of the span-wise direction is indicated,
Figure FDA0003441681810000027
represents the temperature; phi '═ of (ρ', u ', v', w ', T')TFor the disturbance amount, ρ ' represents density, u ' represents flow direction velocity, v ' represents normal velocity, w ' represents span direction velocity, and T ' represents temperature; instantaneous quantity phi and basic flow
Figure FDA0003441681810000028
The equation of the basic flow is subtracted from the equation of the instantaneous quantity, and the nonlinear high-order small quantity is omitted, so that the linear disturbance equation under a Cartesian coordinate system can be obtained, and the form is as follows:
Figure FDA0003441681810000029
wherein the variable φ 'represents the perturbation vector (ρ', u ', v', w ', T')T,Γ、A、B、C、D、Vxx、Vyy、Vzz、Vxy、Vyz、VxzCoefficient matrix representing linear disturbance equation, in which basic flow is contained
Figure FDA00034416818100000210
The information of (a);
3) the linear disturbance equation (2) obtains a linear disturbance equation under a curve coordinate system through coordinate transformation:
Figure FDA00034416818100000211
where ξ, η, ζ denote the flow direction, normal direction and spread direction respectively in a curved coordinate system, and Φ 'denotes a disturbance vector Φ' ((ρ ', u', v ', w', T))T,Γ、
Figure FDA00034416818100000212
D、Vξξ、Vηη
Figure FDA00034416818100000213
Vξη
Figure FDA00034416818100000214
Coefficient matrix representing linear disturbance equation in curve coordinate system, including basic flow
Figure FDA00034416818100000215
The expression of each coefficient matrix is
Figure FDA0003441681810000031
Figure FDA0003441681810000032
Figure FDA0003441681810000033
Vξξ=Vxx,
Figure FDA0003441681810000034
Figure FDA0003441681810000035
Vξη=ηyVxyzVxz,
Vξζ=ζyVzyzVxz,
Vηζ=2ηyζyVyy+2ηzζzVzz+(ηyζzzζy)Vyz
(4)
Wherein J represents Jacobi matrix determinant with expression of
Figure FDA0003441681810000036
4) Deriving a control equation for the matrix-free two-dimensional global stability method, i.e. a frequency-free two-dimensional global stability equation, using a local parallel flow assumption to set the disturbance in the form of a flow direction waveform
Figure FDA0003441681810000037
Wherein
Figure FDA0003441681810000038
A function of the shape representing the time t,
Figure FDA0003441681810000039
wherein
Figure FDA00034416818100000310
The density is expressed as a function of time,
Figure FDA00034416818100000311
the speed of the flow direction is indicated,
Figure FDA00034416818100000312
which is indicative of the normal velocity of the wire,
Figure FDA00034416818100000313
the speed of the span-wise direction is indicated,
Figure FDA00034416818100000314
denotes temperature, α denotes flow direction wave number, c.c. denotes complex conjugation; substituting the formula (6) into a linear disturbance equation (3) under a curve coordinate system, and sorting to obtain a control equation used by the matrix-free global stability method, namely a frequency-free two-dimensional global stability equation:
Figure FDA00034416818100000315
wherein gamma is,
Figure FDA0003441681810000041
Coefficient matrix representing control equation, in which elementary streams are included
Figure FDA0003441681810000042
And a real number form flow direction wave number alpha in a time mode, the expression of each coefficient matrix being
Figure FDA0003441681810000043
2. The method for predicting transition of hypersonic aircraft according to claim 1, wherein in the second step, for equation (7), the relation between the eigenvalue λ and the frequency ω is
Figure FDA0003441681810000044
3. The method for predicting transition of the hypersonic aircraft based on global stability analysis according to claim 1, wherein the determination method of the transition criterion Nt in the step four is as follows: the transition criterion Nt is obtained and corrected through flight experiments or ground wind tunnel experiment data, and the amplitude of disturbance with a certain frequency reaches e of the initial amplitude in the process of propagating towards the flow directionNtThe laminar flow develops into turbulent flow; n obtained in step threemax_allFinding the flow direction position corresponding to Nt on the curve, namely considering the transition starting position Xt of the hypersonic aircraft surface.
CN202110507168.1A 2021-05-10 2021-05-10 Hypersonic aircraft transition prediction method based on global stability analysis Active CN113221350B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110507168.1A CN113221350B (en) 2021-05-10 2021-05-10 Hypersonic aircraft transition prediction method based on global stability analysis

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110507168.1A CN113221350B (en) 2021-05-10 2021-05-10 Hypersonic aircraft transition prediction method based on global stability analysis

Publications (2)

Publication Number Publication Date
CN113221350A CN113221350A (en) 2021-08-06
CN113221350B true CN113221350B (en) 2022-02-18

Family

ID=77094358

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110507168.1A Active CN113221350B (en) 2021-05-10 2021-05-10 Hypersonic aircraft transition prediction method based on global stability analysis

Country Status (1)

Country Link
CN (1) CN113221350B (en)

Families Citing this family (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113657051B (en) * 2021-08-19 2022-05-20 天津大学 Method for predicting surface flow area wall pulsating pressure frequency spectrum of underwater vehicle
CN113704895B (en) * 2021-10-22 2022-10-18 中国空气动力研究与发展中心计算空气动力研究所 Device and method for inhibiting quadratic instability of baby whirlpool and storage medium
CN113998145B (en) * 2022-01-04 2022-03-15 中国空气动力研究与发展中心计算空气动力研究所 Method, device, equipment and medium for detecting instability characteristics of aircraft boundary layer
CN114169267B (en) * 2022-02-11 2022-04-19 中国空气动力研究与发展中心计算空气动力研究所 Method for quickly searching entropy layer characteristic value
CN115146383B (en) * 2022-07-01 2023-01-24 天津大学 Method for forecasting transition position of curved surface boundary layer of super-hydrophobic surface
CN115659522B (en) * 2022-12-27 2023-03-28 中国空气动力研究与发展中心计算空气动力研究所 Aircraft transition position prediction method, device, equipment and medium
CN115859482B (en) * 2023-02-16 2023-04-28 中国空气动力研究与发展中心计算空气动力研究所 Rapid calculation method for steady-state elementary streams of flow field of aircraft
CN116451356B (en) * 2023-05-23 2023-11-24 西安交通大学 Uncertainty compatible natural laminar wing configuration gradient optimization design method
CN118036184B (en) * 2024-02-29 2024-09-17 天津大学 Transition delay optimization method of underwater revolving body head shape based on genetic algorithm
CN118194618B (en) * 2024-05-20 2024-07-16 中国空气动力研究与发展中心计算空气动力研究所 Transition position prediction method
CN118627202B (en) * 2024-07-25 2024-10-11 中国空气动力研究与发展中心计算空气动力研究所 Aircraft transition array plane prediction method based on global stability analysis
CN118656570B (en) * 2024-08-20 2024-10-15 中国空气动力研究与发展中心计算空气动力研究所 Self-similar high-speed boundary layer transition prediction method, device, equipment and medium

Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108287054A (en) * 2017-12-25 2018-07-17 中国航天空气动力技术研究院 A kind of transition Reynolds number acquisition methods under flying condition
CN108304601A (en) * 2017-08-09 2018-07-20 北京空天技术研究所 A kind of judgment method of hypersonic aircraft boundary layer transition
CN108304600A (en) * 2017-08-09 2018-07-20 北京空天技术研究所 A kind of hypersonic aircraft turns to twist position predicting method
CN108304597A (en) * 2017-08-08 2018-07-20 北京空天技术研究所 A kind of high-speed aircraft head leading edge transition prediction device and method
CN108549785A (en) * 2018-05-03 2018-09-18 中国人民解放军国防科技大学 Method for quickly predicting accurate trajectory of hypersonic aircraft based on three-dimensional flight profile
CN110806300A (en) * 2019-10-12 2020-02-18 北京临近空间飞行器系统工程研究所 Measuring point arrangement method suitable for hypersonic flight test transition research
CN111380663A (en) * 2020-02-25 2020-07-07 空气动力学国家重点实验室 Stability method-based cross flow transition experimental data expansion technology
CN111780948A (en) * 2020-06-10 2020-10-16 北京临近空间飞行器系统工程研究所 Method for measuring transition process characteristic of aircraft boundary layer in hypersonic flight test
CN111832159A (en) * 2020-06-23 2020-10-27 北京临近空间飞行器系统工程研究所 Flight test data-based boundary layer transition array surface dynamic evolution process determination method

Family Cites Families (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108216685A (en) * 2016-12-19 2018-06-29 北京空间技术研制试验中心 Suitable for the pneumatic thermal measurement method of blunt body reentry vehicle
CN109033525B (en) * 2018-06-27 2022-08-30 浙江大学 Hypersonic transition prediction method based on simplified three-equation transition model
CN110702356B (en) * 2019-10-12 2021-06-08 空气动力学国家重点实验室 Hypersonic velocity transition prediction method considering surface roughness effect
CN112231847B (en) * 2020-11-04 2024-04-02 中国商用飞机有限责任公司北京民用飞机技术研究中心 Transition position determining method and device, electronic equipment and storage medium

Patent Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108304597A (en) * 2017-08-08 2018-07-20 北京空天技术研究所 A kind of high-speed aircraft head leading edge transition prediction device and method
CN108304601A (en) * 2017-08-09 2018-07-20 北京空天技术研究所 A kind of judgment method of hypersonic aircraft boundary layer transition
CN108304600A (en) * 2017-08-09 2018-07-20 北京空天技术研究所 A kind of hypersonic aircraft turns to twist position predicting method
CN108287054A (en) * 2017-12-25 2018-07-17 中国航天空气动力技术研究院 A kind of transition Reynolds number acquisition methods under flying condition
CN108549785A (en) * 2018-05-03 2018-09-18 中国人民解放军国防科技大学 Method for quickly predicting accurate trajectory of hypersonic aircraft based on three-dimensional flight profile
CN110806300A (en) * 2019-10-12 2020-02-18 北京临近空间飞行器系统工程研究所 Measuring point arrangement method suitable for hypersonic flight test transition research
CN111380663A (en) * 2020-02-25 2020-07-07 空气动力学国家重点实验室 Stability method-based cross flow transition experimental data expansion technology
CN111780948A (en) * 2020-06-10 2020-10-16 北京临近空间飞行器系统工程研究所 Method for measuring transition process characteristic of aircraft boundary layer in hypersonic flight test
CN111832159A (en) * 2020-06-23 2020-10-27 北京临近空间飞行器系统工程研究所 Flight test data-based boundary layer transition array surface dynamic evolution process determination method

Non-Patent Citations (6)

* Cited by examiner, † Cited by third party
Title
An investigation on the disturbance evolution and the transition by resonant-triad interactions with a side-frequency disturbance in a boundary layer;Yongming Zhang etal;《Physics of Fluids》;20200801;第1-25页 *
Application of Arnoldi method to boundary layer instability;Zhang Yong-Ming etal;《Chin. Phys. B》;20151231;第1-13页 *
Extensible Rapid Transition Prediction for Aircraft Conceptual Design;Dev Rajnarayan etal;《29th AIAA Applied Aerodynamics Conference》;20111231;第1-15页 *
一种基于动模态分解的翼型流动转捩预测新方法;韩忠华 等;《航空学报》;20170125;第38卷(第1期);第1-17页 *
高超声速三维边界层横流转捩的数值研究;韩宇峰 等;《空气动力学学报》;20190831;第37卷(第4期);第522-529页 *
高超声速椭圆锥短轴流向涡的二维全局稳定性分析;李晓虎 等;《空气动力学学报》;20180430;第36卷(第2期);第265-272页 *

Also Published As

Publication number Publication date
CN113221350A (en) 2021-08-06

Similar Documents

Publication Publication Date Title
CN113221350B (en) Hypersonic aircraft transition prediction method based on global stability analysis
CN108304600B (en) Method for predicting transition position of hypersonic aircraft
CN108304601B (en) Method for judging transition of boundary layer of hypersonic aircraft
CN111580535B (en) Reentry track three-dimensional profile planning method and system based on convex optimization
Kuzenov et al. Approximate calculation of convective heat transfer near hypersonic aircraft surface
Ding et al. Numerical investigation on the performances of porous matrix with transpiration and film cooling
CN114528665B (en) Prediction method for transition position of underwater flat plate boundary layer heated/cooled on wall surface
CN112861447B (en) Stability theory-based optimal design method for head line type of underwater revolving body
CN113657051B (en) Method for predicting surface flow area wall pulsating pressure frequency spectrum of underwater vehicle
Rodríguez et al. Wavepacket models for subsonic twin jets using 3D parabolized stability equations
Yang et al. An implicit scheme with memory reduction technique for steady state solutions of DVBE in all flow regimes
Kumar et al. Central upwind scheme based immersed boundary method for compressible flows around complex geometries
Kanamori et al. Extension of multipole analysis to laterally asymmetric flowfield around supersonic flight vehicle
Kurbatskii et al. Application of pressure-based coupled solver to the problem of hypersonic missiles with aerospikes
Morimoto et al. Reduction of aerodynamic heating with opposing jet through extended nozzle in high enthalpy flow
Afendikov et al. Numerical simulation of the recirculation flow during the supersonic separation of moving bodies
Liu et al. Nonlinearity analysis of increase in lift of double swept waverider
Chuanzhen et al. Design and analysis of double-swept waverider with wing dihedral
Bedarev et al. Computation of wave interference and relaxation of particles after passing of a shock wave
Bondarev et al. Comparative Estimation of QGDFoam Solver Accuracy for Inviscid Flow Around a Cone
Brahmachary et al. On maximum ballistic coefficient axisymmetric geometries in hypersonic flows
Kim et al. Numerical investigation of a single intermediate-sized bubble in horizontal turbulent channel flow
Kumar et al. Effect of mass injection on secondary instability of hypersonic boundary layer over a blunt cone
Hamilton et al. Improved approximate method for computing convective heating on hypersonic vehicles using unstructured grids
Maochang et al. Numerical Simulation of Gas‐Particle Two‐Phase Flow in a Nozzle with DG 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