Hypersonic transition and forced motion coupling dynamic derivative prediction method

Abstract

The invention provides a hypersonic transition and forced motion coupling dynamic derivative prediction method, and belongs to the technical field of aerospace data processing, and the method comprises the steps: constructing an initial transition model based on a standard Reynolds shear stress transport turbulence model framework by using an intermittent factor transport equation and a transition initial momentum thickness Reynolds number transport equation; carrying out primary correction on a transportation equation of the initial transition model, and carrying out compressibility and transverse flow effect correction on the initial transition model according to a compressible boundary layer theory and transverse flow correlation to obtain a corrected transition model; carrying out coupling calculation on the corrected transition model and forced vibration by utilizing a rigid dynamic grid technology to obtain unsteady aerodynamic force and torque; and identifying the unsteady aerodynamic force and torque by using an integral method to obtain a predicted combined dynamic derivative. The method has a full-dimension dynamic derivative analysis capability, breaks through speed domain limitation, and realizes full coverage from subsonic velocity or transonic velocity or supersonic velocity to hypersonic velocity working conditions.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of aerospace data processing, and particularly relates to a hypersonic transition and forced motion coupled dynamic derivative prediction method. BACKGROUND

[0002] Boundary layer transition refers to the process of fluid in the boundary layer from laminar state to turbulent state, which involves complex fluid dynamics mechanism; the process is jointly affected by multiple factors, and can have significant influence on flow characteristics, aerodynamic force and moment distribution, heat transfer characteristics, etc. Aircraft dynamic characteristics refer to the law of motion state of the aircraft changing with time when the aircraft is subjected to external disturbance or control input. Dynamic derivative is an aerodynamic parameter for describing the dynamic characteristics of the aircraft, indicating the rate of change of aerodynamic force and moment of the aircraft in unsteady motion with respect to motion parameters (such as angular velocity and angular acceleration), and is an important parameter for evaluating the dynamic stability of the aircraft. As a core parameter for the design of aerodynamic dynamic characteristics of the aircraft, dynamic derivative plays a key role in trajectory planning and motion stability evaluation; the parameter is not only a basic aerodynamic index for evaluating the dynamic performance of the aircraft, but also a original basis for constructing the control system, and its accurate prediction is directly related to the reliability of the design and performance evaluation of the aircraft. Before the 1990s, the prediction of dynamic derivative mainly relied on experience and semi-empirical methods, such as modified Newton theory, internal Newton flow theory, Newton-Boltzmann theory and modified shock-expansion wave theory, etc. After the 1990s, with the rapid development of computational fluid dynamics (CFD) technology and computer hardware level, the method of obtaining dynamic derivative based on numerical solution of Euler or Navier-Stokes equation gradually became the mainstream. Domestic and foreign scholars have conducted in-depth research on the method of obtaining dynamic derivative by numerical simulation of forced oscillation or free oscillation process.

[0003] Previous flight test or wind tunnel test researches have shown that boundary layer transition can significantly affect the dynamic characteristics of the aircraft, and produce pitch instability and other phenomena. In the 1960s, researchers observed that when the aircraft has an angle of attack, the transition positions of the leeward surface and the windward surface show asymmetric distribution, and this asymmetry is closely related to the size of the angle of attack, which is called asymmetric transition. In the 1970s, with the in-depth study of slender reentry bodies, it was found that the boundary layer transition flow of the reentry body would cause significant abnormalities in the motion of the cone. This discovery prompted the academic community to focus on the asymmetric transition phenomenon and its influence on the aerodynamic characteristics of the aircraft.

[0004] The research on the influence of transition on the dynamic characteristics of the aircraft started early in China, and existing researchers have The blunt cone model is taken as the research object, a single degree of freedom aerodynamic floating shaft free vibration experiment is carried out in a Ma=6 wind tunnel, and the coupling effect between the unsteady aerodynamic force before and after the boundary layer transition zone moves and the aircraft movement when the model rotates around the pitch axis is mainly studied. The results show that the coupling effect of dynamic transition and pitch motion may cause the pitch stability of the blunt cone aircraft to decrease. The existing researchers take the cone angle The sharp cone Ma=6 flow is taken as the research object, the influence of transition on the static / dynamic stability of the sharp cone under different Reynolds numbers and attack angles is studied through forced vibration dynamic derivative measurement experiment. The results show that the asymmetric transition is the fundamental reason for the change of the static or dynamic stability of the sharp cone, and is closely related to the position of the transition region relative to the center of mass. At present, there is no numerical simulation method for calculating the dynamic derivative of the transition model and forced vibration coupling aircraft. SUMMARY

[0005] In view of the above problems in the prior art, the high supersonic transition and forced motion coupling dynamic derivative prediction method provided by the present application solves the problems of low accuracy of unsteady aerodynamic load calculation, poor dynamic derivative analysis ability, speed domain limitation and difficulty in realizing full coverage of subsonic or transonic or supersonic to hypersonic working conditions in the numerical simulation method of the dynamic derivative of the existing aircraft.

[0006] In order to achieve the above purpose, the technical scheme adopted by the present application is as follows: a high supersonic transition and forced motion coupling dynamic derivative prediction method, comprising the following steps: S1, based on the standard Reynolds shear stress transport turbulence model framework, using the intermittent factor transport equation and the transition starting momentum thickness Reynolds number transport equation, an initial transition model is constructed; S2, the transport equation of the initial transition model is firstly modified, and the compressibility and crossflow effect of the initial transition model are modified according to the compressible boundary layer theory and the crossflow correlation, to obtain the modified transition model; S3, the modified transition model and forced vibration are coupled and calculated by using the rigid dynamic mesh technology, to obtain the unsteady aerodynamic force and moment; S4, based on the pitch forced motion, the unsteady aerodynamic force and moment are identified by using the integral method, to obtain the predicted combined dynamic derivative.

[0007] The beneficial effects of the present application are as follows: the present application establishes an initial transition model, so that the improved transition model is suitable for a high supersonic boundary layer transition model, the dynamic mesh technology is used to realize the coupling calculation of the transition model and forced vibration, the model produces multiple or single degree of freedom steady vibration and maintains the preset frequency and amplitude parameters, and the parameter identification algorithm is combined to obtain the predicted combined dynamic derivative; the dynamic derivative prediction method has full dimensional dynamic derivative analysis ability, breaks through the speed domain limitation, and realizes full coverage of subsonic or transonic or supersonic to hypersonic working conditions; Small amplitude pitching oscillation is described by rigid dynamic mesh technique, and the coupling calculation of the modified transition model and forced oscillation is realized. The relative position between the nodes of the mesh is always kept unchanged in the whole movement process by using the rigid dynamic mesh technique, so that the quality of the mesh is effectively guaranteed without generating new mesh. The calculation amount is greatly reduced without calculating the additional physical quantities such as mesh deformation.

[0008] Further, the S1 comprises the following steps: S101, based on the standard Reynolds shear stress transport turbulence model framework, and according to the experimental data fitting, the intermittency factor transport equation and the transition starting position momentum thickness Reynolds number transport equation are obtained; S102, the separation intermittency factor is obtained by using the transition starting momentum thickness Reynolds number transport equation, and the effective intermittency factor is obtained by using the intermittency factor transport equation; S103, the effective intermittency factor is used to modify the turbulent kinetic energy transport equation in the standard Reynolds shear stress transport turbulence model, and the initial transition model is obtained.

[0009] Further, the expression of the initial transition model is as follows: ; ; ; Wherein, denotes the density, denotes the turbulent fluctuation kinetic energy, denotes the time variable, denotes the velocity component, denotes the spatial variable, denotes the molecular viscosity coefficient, denotes the model constant, denotes the turbulent eddy viscosity coefficient, denotes the generation term of the modified turbulent kinetic energy transport equation, denotes the dissipation term of the modified turbulent kinetic energy transport equation, denotes the effective intermittency factor, denotes the generation term of the turbulent kinetic energy transport equation in the original standard Reynolds shear stress transport turbulence model, denotes the dissipation term of the turbulent kinetic energy transport equation in the original standard Reynolds shear stress transport turbulence model, denotes the intermittency factor, denotes the separation intermittency factor.

[0010] The beneficial effects of the above-mentioned further solutions are as follows: By adding the intermittent factor transport equation and the momentum-thickness-Reynolds number transport equation at the beginning of the transition to the standard Reynolds shear stress transport turbulence model framework, the present invention obtains an initial transition model, making the numerical simulation data of the transition more accurate and having a higher degree of fit with the actual model.

[0011] Furthermore, S2 includes the following steps: S201. By controlling the length of the transition zone, the transport equation of the initial transition model is initially modified, and the local variables are calculated using engineering methods. S202. Based on the compressible boundary layer theory and local variables, a compressible boundary layer correlation function is constructed to analyze the relationship between the pressure gradient parameter and the Mach number, obtain the pressure gradient parameter correction, and introduce the compressible correction into the transition initial momentum thickness Reynolds number transport equation to obtain the compressible correction transition model. S203. Based on the crossflow correlation and local variables, construct the crossflow correlation function, and by modifying the model constants and local variables in the crossflow correlation function, perform crossflow effect correction on the compressible correction transition model to obtain the corrected transition model.

[0012] Furthermore, in step S201, the length of the transition region is controlled by a transition region length function, the expression of which is shown below: ; ; in, The function representing the length of the transition region. This represents the viscous sublayer damping function. Represents the wall vortex Reynolds number. Indicates density, Indicates the distance along the normal direction of the wall. Indicates the turbulent specific dissipation rate. This represents the molecular viscosity coefficient.

[0013] Furthermore, the expression for the compressible boundary layer correlation function is as follows: ; ; in, This represents the compressible boundary layer correlation function. Represents the vortex Reynolds number. This represents the critical momentum thickness Reynolds number. Represents a local function. and Both represent local variables. Indicates the wall temperature. 、 、 All represent fixed parameters.

[0014] Further, the expression of the cross-flow correlation function is as follows: ; Wherein, represents the cross-flow correlation function, represents the cross-flow intensity term, represents the model constant.

[0015] The beneficial effects of the above further scheme are: the initial transition model is modified to control the transition zone length function, and the local variable calculated by the engineering method is used to transform the initial transition model into a new transition model suitable for hypersonic transition prediction, which improves the accuracy of the dynamic derivative prediction. Based on the compressible boundary layer theory, the original pressure gradient parameter is modified by analyzing the relationship between the pressure gradient parameter and the Mach number to obtain the modified pressure gradient parameter, and the compressibility correction is added to the calculation of the original transition starting momentum thickness Reynolds number, realizing the compressibility correction of the transition model. By constructing the cross-flow correlation function and modifying the model constant and the local variable, the cross-flow effect correction of the transition model is realized, and the prediction accuracy of the transition model for hypersonic transition is improved.

[0016] Further, the S4 comprises the following steps: S401, based on the forced motion around the center of mass and the body coordinate system, setting the forced motion form around the center of mass, and obtaining the first and second derivatives of the pitch angle; S402, according to the preset classical dynamic aerodynamic force model, based on the angle of attack of the forced motion form, Taylor expansion is carried out to obtain the Taylor expanded aerodynamic force model; S403, the first and second derivatives of the pitch angle are brought into the Taylor expanded aerodynamic force model, and the linear term is retained to obtain the prediction aerodynamic force model, and the coefficient expression is recorded; S404, according to the coefficient expression, using integral method, by integrating the prediction aerodynamic force model in a cycle, the unsteady aerodynamic force and moment are identified to obtain the prediction combined dynamic derivative.

[0017] Further, the expression of the prediction aerodynamic force model is as follows: ; Wherein, represents the prediction aerodynamic force model, represents the preset classical dynamic aerodynamic force model, represents the Taylor expanded aerodynamic force model, denotes the pitch angle, denotes the first derivative of the pitch angle, denotes the second derivative of the pitch angle, denotes the reduced frequency, denotes the pitch oscillation amplitude, denotes the time variable.

[0018] Further, the predicted combination dynamic derivative expression is as follows: ; ; wherein, denotes the combination static moment dynamic derivative, denotes the sine-weighted periodic integral of the pitch moment, denotes the combination damping dynamic derivative, denotes the cosine-weighted periodic integral of the pitch moment, denotes the initial time of integration, denotes the periodic time.

[0019] The above further scheme has the beneficial effect that the present application uses integral method for parameter identification algorithm, dynamically identifies the predicted combination dynamic derivative, reduces the error probability of prediction, uses periodic integral, the integral term has linear relationship with the oscillation amplitude, the integral is constructed by using the orthogonality of trigonometric functions, separates the aerodynamic force signal and high-frequency noise in the frequency domain, and eliminates the initial state dependence through periodic integral, thereby improving the calculation efficiency and the prediction accuracy. BRIEF DESCRIPTION OF DRAWINGS

[0020] Figure 1 is a flow chart of the method of the present application.

[0021] Figure 2 is a model simulation diagram of the blunt cone calibration model in the present embodiment.

[0022] Figure 3 is a schematic diagram of the geometric shape and the calculation grid of the blunt cone calibration model in the present embodiment.

[0023] Figure 4 is a flow field structure diagram of the blunt cone transition at different times in the present embodiment.

[0024] Figure 5 is a comparison diagram of the calculation results and the experimental results in the present embodiment.

[0025] Figure 6 is a curve diagram of the angle of attack and the moment in the present embodiment.

[0026] Figure 7 is a hysteresis loop diagram of the laminar flow and the transition pitch moment coefficient and the angle of attack in the present embodiment. DETAILED DESCRIPTION

[0027] The specific embodiments of the present application are described below to facilitate the understanding of the present application for those skilled in the art, but it should be clear that the present application is not limited to the scope of the specific embodiments, and for those skilled in the art, it is obvious that various changes are within the spirit and scope of the present application defined and determined by the appended claims, and all the inventions utilizing the concept of the present application are within the scope of protection.

[0028] Before explaining the present embodiment, the following terms are explained: SST turbulent flow model framework: standard Reynolds shear stress transport turbulent flow model framework; Etkin aerodynamic force model: a model for describing the aerodynamic force hysteresis effect in the dynamic motion of the aircraft; Blasius boundary layer: incompressible boundary layer on a flat plate without pressure gradient; CFD: computational fluid dynamics.

[0029] Embodiment 1 As shown in Figure 1 The present application provides a hypersonic transition and forced motion coupled dynamic derivative prediction method, which is realized as follows: S1, based on the standard Reynolds shear stress transport turbulent flow model framework, using the intermittency factor transport equation and the transition starting momentum thickness Reynolds number transport equation, an initial transition model is constructed, the specific steps are as follows: S101, based on the standard Reynolds shear stress transport turbulent flow model framework, and according to the experimental data fitting, the intermittency factor transport equation and the momentum thickness Reynolds number transport equation at the transition starting position are obtained; S102, using the transition starting momentum thickness Reynolds number transport equation, the separation intermittency factor is obtained, and using the intermittency factor transport equation, the effective intermittency factor is obtained; S103, using the effective intermittency factor, the turbulent kinetic energy transport equation in the standard Reynolds shear stress transport turbulent flow model is modified, and the initial transition model is obtained.

[0030] In this embodiment, the initial transition model is γ Reθ transition model, based on the SST turbulent flow model framework, two additional transport equations, intermittency factor transport equation and momentum thickness Reynolds number transport equation at the transition starting position, are added, the above two transport equations are based on experimental data fitting, the specific expression is as follows: ; ; Where, ρ represents the density, η represents the intermittency factor, denotes a time variable, denotes a velocity component, denotes a spatial variable, denotes a molecular viscosity coefficient, denotes a turbulent eddy viscosity coefficient, and denote model constants, denotes a generation term of the intermittency factor transport equation, denotes a dissipation term of the intermittency factor transport equation, denotes a momentum thickness Reynolds number at the transition onset location, denotes a generation term of the momentum thickness Reynolds number transport equation at the transition onset location; the generation term of the intermittency factor transport equation and the dissipation term are expressed as follows: ; ; wherein, denotes a transition region length function for controlling the transition region length, denotes a modulus of the strain rate tensor, denotes a switch function triggering the transition, denotes a vorticity tensor, denotes a dissipation function for controlling the intermittency factor ; , , , denote model constants; the generation term of the momentum thickness Reynolds number transport equation at the transition onset location is expressed as follows: ; wherein, denotes a fixed parameter, being 0.03, denotes a new transition momentum thickness Reynolds number in the free stream, denotes a boundary layer indicator, being 1 when inside the boundary layer and 0 when outside the boundary layer, the parameter being used to switch off the source term of the transport equation, only allowing the diffusion process in the boundary layer.

[0031] In the embodiment, the separation intermittency factor, γ is obtained by using the transition onset momentum thickness Reynolds number transport equation, and the effective intermittency factor is obtained by using the intermittency factor transport equation; the expression of the separation intermittency factor is as follows: ; wherein, is a constant that can be adjusted by tuning the size of to adjust the sensitivity of the model to separation-induced transition, is the vortex Reynolds number, is the critical momentum thickness Reynolds number, is the separation reattachment function; The expression of the effective intermittency factor is as follows: ; The generation term and the dissipation term of the turbulent kinetic energy transport equation in the original SST turbulence model are modified using the effective intermittency factor to realize the coupling with the SST turbulence model, and an initial transition model is obtained; the expression of the initial transition model is as follows: ; ; wherein, is the turbulent fluctuation kinetic energy, is the modified generation term of the turbulent kinetic energy transport equation, is the modified dissipation term of the turbulent kinetic energy transport equation, is a model constant; The transport equation expression of the specific dissipation rate in the original SST turbulence model is as follows: ; wherein, is the turbulent specific dissipation rate, is the generation term of the specific dissipation rate transport equation, is the dissipation term of the specific dissipation rate transport equation, is a mixing function, which is 1 in the boundary layer and 0 in the free stream, and both represent model constants.

[0032] S2, the transport equation of the initial transition model is modified for the first time, and the compressibility and crossflow effect of the initial transition model are modified according to the compressible boundary layer theory and the crossflow correlation, to obtain a modified transition model, and the specific steps are as follows: S201, the transport equation of the initial transition model is modified for the first time by controlling the length of the transition zone, and the local variables are calculated by using an engineering method. S202, according to the compressible boundary layer theory and local variables, the compressible boundary layer correlation function is constructed, the relationship between the pressure gradient parameter and the Mach number is analyzed, the pressure gradient parameter correction is obtained, and the compressible correction is introduced into the transition starting momentum thickness Reynolds number transport equation, and the transition model corrected by the compressible correction is obtained.

[0033] In this embodiment, the γ Reθtransition model is modified to have a new criterion suitable for hypersonic transition prediction, and the initial modified γ Reθtransition model is obtained by controlling the transition zone length. The transition zone length function controls the transition zone length, and the expression is as follows: ; ; Wherein, The viscous sublayer damping function is represented by The wall vortex Reynolds number is represented by The normal distance along the wall is represented by And the local variables And The local variable Determines the non-local Mach number Of the boundary layer edge.

[0034] In this embodiment, the compressibility of the initial modified γ Reθtransition model is corrected, and the high-speed boundary layer local correlation function based on the compressible boundary layer self-similar solution is constructed according to the compressible boundary layer theory and local variables, and the specific expression is as follows: ; ; Wherein, The compressible boundary layer correlation function is represented by The vortex Reynolds number is represented by The critical momentum thickness Reynolds number is represented by The local function is represented by And Both represent local variables, The wall temperature is represented by , , Both represent parameters calculated by local variables, and the calculation expression is as follows:

[0035] Based on the compressible boundary layer theory, by analyzing the relationship between the pressure gradient parameter and the Mach number, the following expression of the pressure gradient parameter correction is obtained: ; wherein, represents the corrected pressure gradient parameter, represents the original pressure gradient parameter, represents the specific heat ratio, and in the embodiment, the specific heat ratio is 1.4; In the calculation of the original transition starting momentum thickness Reynolds number , the compressibility correction is added, and by obtaining the corrected original transition starting momentum thickness Reynolds number , the transition model corrected by compressibility is obtained. The expression of the compressibility correction is as follows: ; ; wherein, represents the corrected original transition starting momentum thickness Reynolds number, represents the original transition starting momentum thickness Reynolds number.

[0036] S203, according to the cross-flow correlation and the local variable, a cross-flow correlation function is constructed, and by correcting the model constant in the cross-flow correlation function and the local variable, the transition model corrected by compressibility is corrected by the cross-flow effect, and the corrected transition model is obtained.

[0037] In the embodiment, the cross-flow correlation is based on the concept of local helicity, which is widely used as a cross-flow indicator in the model based on local correlation, and a cross-flow correlation function is constructed, and the expression is as follows: ; wherein, represents the cross-flow correlation function, represents the cross-flow intensity term, represents the model constant, which is calibrated to 28.0 in a noise environment, and is calibrated to 45.0 in a quiet environment (the free flow turbulence intensity is much smaller than 0.1%); the cross-flow intensity term expression is constructed as follows: ; ; ; wherein, represents the dimensionless cross-flow intensity, which can be defined by the local helicity, represents the local velocity magnitude, represents the distance along the wall normal direction, denotes local helicity; The transition model is modified by modifying the model constant and local variable in the cross-flow correlation function, and the transition model is modified by modifying the cross-flow effect. The starting function of the modified transition model is defined as follows: ; Wherein, denotes the primary transition trigger function, denotes the fourth power enhancement term of the primary trigger function, denotes the enhanced transition criterion, denotes the turbulence shielding function, denotes the switch function of triggering transition, which is the final transition starting function.

[0038] S3, the modified transition model and forced vibration are coupled by using rigid dynamic grid technology, and the unsteady aerodynamic force and torque are obtained.

[0039] In this embodiment, the rigid dynamic grid technology is used to describe the small amplitude pitch vibration, and the coupling calculation of the modified transition model and the forced vibration is realized to obtain the unsteady aerodynamic force and torque. The idea of rigid dynamic grid technology is to rigidly connect the calculation grid with the physical shape and move together. The unsteady aerodynamic force and torque are calculated by using the solver in computational fluid dynamics. Compared with the non-structural grid reconstruction technology, the grid reconstruction needs a large amount of calculation resources, time and effort. For example, the sliding grid technology needs to build different sliding surfaces when simulating different types of motion, and different grids need to be generated. It is also time-consuming and laborious to calculate large-scale conditions. The rigid dynamic grid technology has obvious advantages: the grid does not need to be regenerated during the whole unsteady calculation process, is simple and suitable for any amplitude of motion, and the grid topology structure does not change during the unsteady motion process, so different grids do not need to be generated for different motion types, realizing the separation of grid and motion. In the whole motion process, the relative position between the grid nodes remains unchanged, so there is no need to generate new grids, which can effectively ensure the quality of the grid. Moreover, since no additional physical quantities such as grid deformation need to be calculated, the amount of calculation is greatly reduced. This is an advantage for the numerical simulation of the coupling of boundary layer transition and forced motion by this method.

[0040] S4, based on the pitch forced motion, the unsteady aerodynamic force and torque are identified by using integral method to obtain the predicted combined dynamic derivative, the specific steps are as follows: S401, based on the pitch forced motion and the body coordinate system, the forced motion form around the center of mass is set, and the first and second derivatives of the pitch angle are obtained; S402, according to the preset classical dynamic aerodynamic force model, based on the angle of attack of forced motion form, Taylor expansion is carried out to obtain the Taylor expanded aerodynamic force model; S403, the first order derivative and the second order derivative of the pitch angle are brought into the Taylor expanded aerodynamic force model, and the linear term is retained to obtain the predicted aerodynamic force model, and the coefficient expression is recorded; S404, according to the coefficient expression, the integral method is used to identify the unsteady aerodynamic force and moment by integrating the predicted aerodynamic force model in a cycle to obtain the predicted combined dynamic derivative.

[0041] In this embodiment, based on the pitch forced motion, the integral method is used to identify the unsteady aerodynamic moment to obtain the predicted combined dynamic derivative; the preset classical dynamic aerodynamic force model is Etkin classical dynamic aerodynamic force model; The pitch damping derivative is the total damping derivative of the pitch channel, which is a combined derivative, including the direct damping derivative And the washout time difference derivative , which directly reflects the dynamic stability of the pitch channel, which is called static derivative, pitch damping derivative; The identification by using the integral method can be carried out by using single point method, integral method and least square method; In the body axis coordinate system, the forced motion form (dimensionless form) around the center of mass is set as follows: ; ; Among them, indicates the power angle, indicates the pitch angle, indicates the reduced frequency, indicates the oscillation angular frequency, indicates the reference length, indicates the flow velocity; And the pitch angle is derived to obtain the first order derivative and the second order derivative of the pitch angle, and the expression is as follows: ; Among them, indicates the first order derivative of the pitch angle, indicates the second order derivative of the pitch angle, indicates the pitch oscillation amplitude; considering the longitudinal motion of the aircraft, if the size of the center of mass velocity is limited to be constant, according to the Etkin classical dynamic aerodynamic force model, the expression is as follows: ; Among them, indicates the pitch moment coefficient at the moment This represents the derivative of the pitch moment coefficient with respect to the angle of attack. Represents a time variable. express The change in angle of attack relative to the equilibrium position at any given moment. express The change in the rate of change of angle of attack relative to the equilibrium position at any given time. express The change in pitch angular velocity relative to the equilibrium position at any given time. express The change in the rate of change of pitch angular velocity relative to the equilibrium position at any given moment; Furthermore, a Taylor expansion is performed on the angle of attack based on the forced motion form, resulting in the Taylor-expanded aerodynamic model, expressed as follows: ; in, , The subscript 0 indicates the angle of attack. The value at that location.

[0042] In this embodiment, the aerodynamic model expanded by Taylor can also be written as: ; in, Defined as angle of attack The pitch static derivative at that point, Defined as the pitch derivative. Defined as the second-order term of the pitch motion derivative, and substituting the expressions for the first and second derivatives of the pitch angle into the aerodynamic model after Taylor expansion, while retaining the linear terms, the predicted aerodynamic model is obtained, and the coefficient expressions are recorded; the expression of the predicted aerodynamic model is as follows: ; in, This indicates a predictive aerodynamic model. This indicates a pre-defined classical dynamic aerodynamic model. This represents the aerodynamic model expanded by Taylor; The expression for the recording coefficient is: ; Furthermore, by using the integral method, the unsteady aerodynamic forces and torques are identified by integrating the predicted aerodynamic model over one period, and the predicted combined dynamic derivatives are obtained. The expression for the predicted combination dynamic derivative is as follows: ; ; in, Represents the dynamic derivative of the combined static moment. represents the sine-weighted periodic integral of the pitch moment, represents the combined damping dynamic derivative, represents the cosine-weighted periodic integral of the pitch moment, represents the initial time of integration, represents the periodic time; represents the preset classical dynamic aerodynamic force model, as a moment coefficient in prediction, after simulation, the solver can directly output specific values, since forced motion is performed at this time, the angle of attack changes between -2~+2, each angle of attack can output a moment coefficient, and can obtain the curve function changing with the angle of attack , the longitudinal combined dynamic derivative is obtained by integrating the unsteady aerodynamic force and moment identification; the aerodynamic force and moment are dimensionless, and the moment coefficient and the aerodynamic force coefficient are taken to participate in prediction.

[0043] In this embodiment, a hypersonic transition and forced motion coupled dynamic derivative prediction method is provided, the transition model adopted is a γ Re transition model improved for hypersonic boundary layer transition model, then the dynamic mesh technology is applied to realize the coupling calculation of the transition model and the forced vibration, so that the model produces multi / single degree of freedom steady vibration and maintains the preset frequency and amplitude parameters. After numerical simulation, the response data of each degree of freedom are obtained, and the dynamic stability derivative can be obtained by combining the parameter identification algorithm. The comparison results of numerical calculation and experiment show that the method has very high accuracy. Corresponding to the forced vibration method in wind tunnel experiment technology, the dynamic parameter identification CFD method becomes a feasible technical path for dynamic derivative calculation. The method realizes the accurate calculation of unsteady aerodynamic load by setting the forced vibration working conditions of parameters such as amplitude and frequency.

[0044] Embodiment 2 In this embodiment, as shown in Figure 2 and Figure 3 , the calculation model selects the HSCM-3 blunt cone standard model specified in GJB-4399-2002 specification, and the geometric characteristic parameters are: half cone angle = 5°, cone length L = 421 mm, base diameter = 148.3 mm, head radius = 8.75 mm. The configuration has typical blunt head characteristics, and the blunt ratio = 0.118, and the mass center of the model is located at the cone length = 64%. The reference cone length is L = 0.421 m, and the reference area is the bottom cross-sectional area S = 0.0173 ​​​​​​​. Calculation conditions: Mach number = 6, forced oscillation amplitude A = 2°, dimensionless reduced frequency = 0.1, where = 0.1, the calculation can be obtained = 37.2 Hz, calculation time step Timestep = 0.0001, inner iteration 50 steps; under the above conditions, dynamic transition and laminar flow dynamic contrast calculation is carried out, and the calculation conditions are shown in Table 1.

[0045] Table 1

[0046] In Table 1, Case represents the example, Re ∞ / m represents the unit Reynolds number, AOA represents the attack angle, T ∞ ,K represents the incoming flow temperature, Tu ∞ ,% represents the flow turbulence.

[0047] In this embodiment, as Figure 4 shown, the flow field Mach number and intermittent factor calculation results of the transition model at different times can be seen. The Mach number cloud map of the symmetry plane at t = 0s, and the wall surface is the intermittent factor distribution. Under the action of high Mach number strong compression effect airflow, the blunt leading edge will form a significant detached bow shock structure. The moments t = 0.079s and t = 0.124s correspond to the flow field at the attack angle and respectively. By comparing, it is found that under the condition of non-zero attack angle, the wall surface streamline is curved, there is obvious crossflow effect, and from the wall surface intermittent factor distribution, it can also be found that the crossflow induces transition.

[0048] In this embodiment, the transverse flow is the key control factor for inducing the transition of the blunt cone boundary layer, as Figure 5 shown, the comparison results of the wall heat flow distribution under laminar and transition states at t = 0.229 s are shown. The numerical calculation results and experimental observation data are highly consistent in the transition front position (deviation < 2%), which effectively verifies the reliability of the dynamic transition model and the forced motion coupling solving method. As Figure 6 shown, the curves of the attack angle and the pitch moment with time reveal significant unsteady dynamic characteristics: there is obvious phase lag between the valley value of the attack angle and the moment peak value calculated by the transition model in the same period, and the phase lag between the moment peak value calculated by the transition model and the moment peak value calculated by the laminar flow is relatively small. The difference mechanism essentially reflects the transient effect in the flow separation and vortex structure evolution process. As Figure 7As shown, the combined pitch moment coefficient hysteresis loops of laminar and transition conditions based on the integral identification method all present anticlockwise rotation characteristics. According to the aeroelastic theory, the rotation direction of the hysteresis loop determines the sign of the dynamic stability derivative, while the enclosed area represents its magnitude. The quantitative comparison shows that the hysteresis loop area of transition condition is reduced by 35.55% compared with that of laminar condition, and its ellipticity parameter decreases from 1.25 to 0.92. This geometric characteristic variation reveals that the boundary layer transition induced by the crossflow will significantly weaken the pitch damping characteristics of the vehicle; The static stability derivative and dynamic stability derivative results identified by the integral method are shown in Table 2.

[0049] Table 2

[0050] It can be seen that the coupling effect of transition induced by crossflow and forced motion in the pitch channel will lead to the decrease of static / dynamic stability of the blunt cone. This is consistent with the conclusion obtained in the experiment, and further proves the correctness of the dynamic transition and forced motion coupling calculation method in this paper.

Claims

1. A method of hypersonic transition and forced motion coupled dynamic derivative prediction, characterized in that, The method comprises the following steps: S1, based on the standard Reynolds shear stress transport turbulence model framework, using the intermittency factor transport equation and the transition starting momentum thickness Reynolds number transport equation, an initial transition model is constructed; S2, the transport equation of the initial transition model is first modified, and according to the compressible boundary layer theory and the crossflow correlation, the compressibility and crossflow effect of the initial transition model are modified to obtain the modified transition model; S3, using the rigid dynamic mesh technology, the modified transition model and forced vibration are coupled to obtain unsteady aerodynamic force and moment; S4, based on the forced pitching motion, the unsteady aerodynamic force and moment are identified by using integral method to obtain the predicted combined dynamic derivative.

2. The coupled hypersonic transition and forced motion dynamic derivative prediction method of claim 1, wherein, The S1 comprises the following steps: S101, based on the standard Reynolds shear stress transport turbulence model framework, and according to the experimental data fitting, the intermittency factor transport equation and the transition starting position momentum thickness Reynolds number transport equation are obtained; S102, using the transition starting momentum thickness Reynolds number transport equation, the separation intermittency factor is obtained, and using the intermittency factor transport equation, the effective intermittency factor is obtained; S103, using the effective intermittency factor, the turbulent kinetic energy transport equation in the standard Reynolds shear stress transport turbulence model is modified to obtain the initial transition model.

3. The coupled hypersonic transition and forced motion dynamic derivative prediction method of claim 2, wherein, The expression of the initial transition model is as follows: wherein denotes the density, denotes the turbulent fluctuation kinetic energy, denotes the time variable, denotes the velocity component, denotes the spatial variable, denotes the molecular viscosity coefficient, denotes the model constant, denotes the turbulent eddy viscosity coefficient, denotes the production term of the modified turbulent kinetic energy transport equation, denotes the dissipation term of the modified turbulent kinetic energy transport equation, denotes the effective intermittency factor, denotes the production term of the turbulent kinetic energy transport equation in the original standard Reynolds shear stress transport turbulence model, denotes the dissipation term of the turbulent kinetic energy transport equation in the original standard Reynolds shear stress transport turbulence model, denotes the intermittency factor, denotes the separation intermittency factor.

4. The coupled hypersonic transition and forced motion dynamic derivative prediction method of claim 1, wherein, The S2 comprises the following steps: S201, the transport equation of the initial transition model is first modified by controlling the transition zone length, and the local variable is calculated by using the engineering method; S202, according to the compressible boundary layer theory and the local variable, the compressible boundary layer correlation function is constructed, the relationship between the pressure gradient parameter and the Mach number is analyzed, the pressure gradient parameter is modified, and the compressible correction is introduced into the transition starting momentum thickness Reynolds number transport equation to obtain the transition model after compressible correction; S203, according to the crossflow correlation and the local variable, the crossflow correlation function is constructed, and the model constant in the crossflow correlation function and the local variable are modified to modify the transition model after compressible correction to obtain the modified transition model.

5. The coupled hypersonic transition and forced motion dynamic derivative prediction method of claim 4, wherein, In the S201, the transition zone length is controlled by the transition zone length function, and the expression of the transition zone length function is as follows: wherein, represents a transition region length function, represents a viscous sublayer damping function, represents a wall vortex Reynolds number, represents a density, represents a distance along the wall normal, represents a turbulent specific dissipation rate, represents a molecular viscosity coefficient.

6. The coupled hypersonic transition and forced motion dynamic derivative prediction method of claim 5, wherein, The expression of the compressible boundary layer correlation function is as follows: wherein, represents the compressible boundary layer correlation function, represents the vortex Reynolds number, represents the critical momentum thickness Reynolds number, represents the local function, and both represent the local variable, represents the wall temperature, , , both represent the fixed parameter.

7. The coupled hypersonic transition and forced motion dynamic derivative prediction method of claim 6, wherein, The expression of the crossflow correlation function is as follows: wherein, represents the cross-flow correlation function, represents the cross-flow intensity term, represents the model constant.

8. The coupled hypersonic transition and forced motion dynamic derivative prediction method of claim 1, wherein, The S4 comprises the following steps: S401, based on the forced pitching motion and the body coordinate system, the forced motion form around the center of mass is set, and the first and second derivatives of the pitch angle are obtained; S402, according to the preset classical dynamic aerodynamic force model, based on the angle of attack of the forced motion form, Taylor expansion is carried out to obtain the Taylor expanded aerodynamic force model; S403, the first and second derivatives of the pitch angle are brought into the Taylor expanded aerodynamic force model, and the linear term is retained to obtain the predicted aerodynamic force model, and the coefficient expression is recorded; S404, according to the coefficient expression, using integral method, by integrating the predicted aerodynamic force model in a period, the unsteady aerodynamic force and moment are identified to obtain the predicted combined dynamic derivative.

9. The coupled hypersonic transition and forced motion dynamic derivative prediction method of claim 8, wherein, The expression of the predicted aerodynamic force model is as follows: wherein, denotes a predicted aerodynamic force model, denotes a pre-set classical dynamic aerodynamic force model, denotes a Taylor expanded aerodynamic force model, denotes a pitch angle, denotes a pitch angle first derivative, denotes a pitch angle second derivative, denotes a reduced frequency, denotes a pitch oscillation amplitude, denotes a time variable.

10. The coupled hypersonic transition and forced motion dynamic derivative prediction method of claim 9, wherein, The expression of the predicted combined dynamic derivative is as follows: wherein, denotes the combined moment derivative, denotes the sine weighted periodic integral of the pitch moment, denotes the combined damping derivative, denotes the cosine weighted periodic integral of the pitch moment, denotes the integration initial time, denotes the periodic time.

Images