Hypersonic flow instability multi-field coupling prediction method and system

Abstract

The invention discloses a hypersonic flow instability multi-field coupling prediction method and system, and relates to the technical field of modern high-speed aircraft design, and the method comprises the steps: obtaining the multi-physical field parameter distribution of a hypersonic flow field; based on a shock wave-boundary layer interaction mechanism, constructing an instability sensitive graph data structure by taking multi-physical field parameters as input, setting graph nodes in a shock wave influence region, a boundary layer region and an intersection region thereof, and calculating connection weights among the nodes; performing multi-scale decomposition on the instability sensitive graph data structure, and extracting instability modes of a shock wave scale, a vortex scale and a molecular scale; screening the instability modes by adopting a hypersonic velocity stability criterion, and determining a dominant instability mode; and constructing an instability evolution prediction model based on the dominant instability mode, and outputting a spatio-temporal evolution rule of hypersonic flow instability. According to the method, the problem that the flow instability phenomenon in the hypersonic velocity complex multi-physical field environment is difficult to accurately predict can be solved.

Description

Technical Field

[0001] The present invention relates to the technical field of modern high-speed aircraft design, in particular to a method and system for predicting multi-field coupling of hypersonic flow instability. Background Art

[0002] Modern high-speed aircraft design technology is a comprehensive field that involves the intersection and integration of multiple disciplines and technologies. It aims to design aircraft that can operate effectively under hypersonic conditions. Therefore, how to use advanced technical means to improve the intelligence level and safety of modern high-speed aircraft design has become one of the urgent issues to be solved.

[0003] In the field of modern high-speed aircraft design, traditional single-temperature models cannot accurately describe the non-equilibrium energy exchange between electrons and heavy particles in plasma under hypersonic conditions, resulting in inaccurate simulation results. Existing fluid mechanics and electromagnetic field analysis methods are usually performed independently, lacking effective multi-physics field coupling analysis methods, making it difficult to fully understand the physical processes in high-speed flow environments. At the same time, traditional numerical simulation methods have difficulty handling data integration and analysis at different scales from micro to macro. Summary of the Invention

[0004] In view of the above existing problems, the present invention is proposed.

[0005] Therefore, the present invention provides a multi-field coupling prediction method for hypersonic flow instability to solve the problem that flow instability phenomena in a complex hypersonic multi-physical field environment are difficult to accurately predict.

[0006] In order to solve the above technical problems, the present invention provides the following technical solutions:

[0007] In a first aspect, the present invention provides a multi-field coupling prediction method for hypersonic flow instability, which comprises:

[0008] S1. Obtaining the multi-physics parameter distribution of the hypersonic flow field;

[0009] S2. Based on the shock wave-boundary layer interaction mechanism, the instability sensitivity map data structure is constructed with multi-physics field parameters as input. Graph nodes are set in the shock wave impact area, boundary layer area, and their intersection area, and the connection weights between nodes are calculated.

[0010] S3. Perform multi-scale decomposition on the instability sensitivity map data structure to extract the instability modes at the shock wave scale, vortex scale, and molecular scale;

[0011] S4. Use hypersonic stability criteria to screen the instability modes and determine the dominant instability mode;

[0012] S5. Construct an instability evolution prediction model based on the dominant instability mode and output the spatiotemporal evolution law of hypersonic flow instability.

[0013] Preferably, the method for setting the graph node in step S2 includes:

[0014] Setting a shock wave instability node in an area where the pressure gradient exceeds a first preset threshold;

[0015] Setting a shear instability node in an area where the velocity gradient exceeds a second preset threshold;

[0016] Set up a coupling instability node in the intersection area of ​​shock wave and boundary layer;

[0017] The instability strength index of each node is calculated according to the following formula:

[0018]

[0019] in, is the pressure gradient modulus, δ is the boundary layer thickness, p avg is the local average pressure, is the normal velocity gradient, U ∞ is the incoming flow velocity, T e is the electron temperature, T h is the heavy particle temperature.

[0020] Preferably, the connection weights between nodes in step S2 are determined based on the physical transmission mechanism, and the weight calculation formula is:

[0021] w ij =W conv ·f conv (Pe ij ,|r i -r j |)+W diff ·f diff (Pe ij ,|r i -r j |)+W em

[0022] ·f em (Re m,ij ,|r i -r j |)

[0023] in, is the Peclet number, which characterizes the competition intensity between convection and diffusion, u i ,u j is the speed of nodes i and j, r i ,r j is the node position, v is the kinematic viscosity; is the magnetic Reynolds number, is the magnetic diffusivity, σ e is the electrical conductivity, μ0 is the vacuum magnetic permeability; W conv 、W diff 、W em are the weight coefficients of convective transmission, diffusion transmission and electromagnetic transmission respectively; f conv 、f diff 、f em are the corresponding convective transfer function, diffusion transfer function, and electromagnetic transfer function, respectively.

[0024] Preferably, the specific form of each transfer function is:

[0025]

[0026]

[0027] is the convection characteristic length, Δt is the time step, and CFL is the dimensionless CFL number; is the diffusion characteristic length; is the characteristic length of plasma, c is the speed of light, ω p is the plasma frequency.

[0028] Preferably, the multi-scale decomposition in step S3 includes:

[0029] Construct Laplace operators at three scale levels: shock-wave scale operators are used to capture pressure wave propagation characteristics, vortex-scale operators are used to capture the evolution characteristics of rotating structures, and molecular-scale operators are used to capture viscous diffusion characteristics.

[0030] Perform eigenvalue decomposition on the Laplace operator of each scale to obtain the unstable mode of the corresponding scale;

[0031] The modes at the three scales are merged to form a complete multi-scale unstable mode set.

[0032] Preferably, the hypersonic stability criterion in step S4 includes a critical growth rate determination, and the critical growth rate is calculated as follows:

[0033]

[0034] Among them, Ma is the Mach number, Re is the Reynolds number, and Ha is the Hartmann number.

[0035] Preferably, the determination of the dominant instability mode in step S4 is based on a comprehensive evaluation of the modal growth rate, spatial influence range and energy contribution.

[0036] Preferably, the instability evolution prediction model in step S5 includes:

[0037] The instability evolution equation based on the linear superposition of the dominant modes is established in the following form:

[0038]

[0039] Among them, k is the time-dependent amplitude of the kth dominant instability mode, λ k is the linear growth rate, N kmn is the nonlinear coupling coefficient;

[0040] The spatiotemporal distribution prediction results of the instability field variables are obtained by solving the evolution equation.

[0041] Preferably, the multi-physical field parameters in step S1 include density field, velocity field, pressure field, temperature field, electron temperature field and magnetic field intensity.

[0042] In a second aspect, the present invention provides a hypersonic flow instability multi-field coupling prediction system, comprising:

[0043] A data acquisition module is used to acquire and pre-process the multi-physics parameters of the hypersonic flow field;

[0044] A graph construction module is used to set graph nodes in the shock wave impact area, boundary layer area, and their intersection area according to preset thresholds, calculate the node instability intensity index, and calculate the connection weights between nodes based on convective transmission, diffusion transmission, and electromagnetic transmission mechanisms to construct the instability sensitivity graph data structure;

[0045] The multi-scale decomposition module is used to construct feature extraction operators at the shock wave scale, vortex scale, and molecular scale, perform eigenvalue decomposition on each scale operator, extract the corresponding unstable modes, and merge them to form a complete unstable mode set;

[0046] Mode screening module, used to determine the dominant instability mode based on critical growth rate determination and comprehensive importance evaluation;

[0047] The prediction module is used to construct an instability evolution prediction model including nonlinear interactions based on the dominant instability mode, and solve the evolution equation to output the spatiotemporal distribution law of instability.

[0048] The beneficial effects of the present invention are:

[0049] (1) The present invention achieves a systematic prediction of hypersonic multi-physics field coupled flow instability phenomena through the organic combination of instability sensitivity map data structure and multi-scale decomposition technology. By setting targeted nodes in the shock wave impact area, boundary layer area and their intersection area and calculating the connection weight based on physical mechanisms, the computational efficiency is improved compared with the traditional full-domain grid method, while ensuring the accurate capture of key instability areas, providing a computationally feasible and physically reasonable prediction method for complex hypersonic flow instability problems.

[0050] (2) By introducing a node instability strength index that integrates pressure gradient, velocity gradient, and the difference between electron temperature and heavy particle temperature, a quantitative characterization of the multi-physics coupling effect is achieved. Combined with the connection weight calculation based on the Peclet number and magnetic Reynolds number, the competing relationships among convective transport, diffusive transport, and electromagnetic transport mechanisms are effectively captured, enabling the constructed instability sensitivity map to accurately reflect the physical transport characteristics and instability sensitivity distribution in the hypersonic flow field.

[0051] (3) The Laplace operator at three scale levels is used to capture the characteristics of pressure wave propagation, rotational structure evolution, and viscous diffusion respectively. The instability modes corresponding to each scale are extracted through eigenvalue decomposition, effectively separating the instability phenomena dominated by different physical mechanisms. This multi-scale processing method avoids the limitations of single-scale analysis and can comprehensively identify the instability mechanisms at all levels from the macroscopic shock wave structure to the microscopic diffusion process, thereby improving the completeness and accuracy of modal identification.

[0052] (4) By introducing a critical growth rate determination formula that takes into account the Mach number, Reynolds number, and Hartmann number, combined with a comprehensive evaluation of modal growth rate, spatial influence range, and energy contribution, the goal of accurately screening the dominant instability mode from a large number of candidate modes is achieved;

[0053] (5) An instability evolution equation containing linear growth terms and nonlinear coupling terms is established. The interaction between different modes is described by the nonlinear coupling coefficient, which can accurately predict the complete evolution process of the instability phenomenon from linear development to nonlinear saturation. BRIEF DESCRIPTION OF THE DRAWINGS

[0054] In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the following briefly introduces the drawings required for use in the description of the embodiments. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0055] Figure 1 is a flow chart of the method of the present invention;

[0056] Figure 2 It is a technical implementation diagram of the present invention;

[0057] Figure 3 This is a system structure diagram of the present invention. DETAILED DESCRIPTION

[0058] In order to make the above-mentioned objects, features and advantages of the present invention more obvious and easy to understand, the specific embodiments of the present invention are described in detail below with reference to the accompanying drawings.

[0059] In the following description, many specific details are set forth to facilitate a full understanding of the present invention. However, the present invention may also be implemented in other ways different from those described herein. Those skilled in the art may make similar generalizations without violating the connotation of the present invention. Therefore, the present invention is not limited to the specific embodiments disclosed below.

[0060] Secondly, the term "one embodiment" or "embodiment" herein refers to a specific feature, structure, or characteristic that may be included in at least one implementation of the present invention. The phrase "in one embodiment" appearing in various places throughout this specification does not necessarily refer to the same embodiment, nor does it refer to a separate or selective embodiment that is mutually exclusive of other embodiments.

[0061] like Figure 1 As shown, the present invention provides a multi-field coupling prediction method for hypersonic flow instability, comprising the following steps:

[0062] S1. Obtaining the multi-physics parameter distribution of the hypersonic flow field;

[0063] S2. Based on the shock wave-boundary layer interaction mechanism, the instability sensitivity map data structure is constructed with multi-physics field parameters as input. Graph nodes are set in the shock wave impact area, boundary layer area, and their intersection area, and the connection weights between nodes are calculated.

[0064] S3. Perform multi-scale decomposition on the instability sensitivity map data structure to extract the instability modes at the shock wave scale, vortex scale, and molecular scale;

[0065] S4. Use hypersonic stability criteria to screen the instability modes and determine the dominant instability mode;

[0066] S5. Construct an instability evolution prediction model based on the dominant instability mode and output the spatiotemporal evolution law of hypersonic flow instability.

[0067] like Figure 2 As shown, the technical implementation process of the present invention is as follows: first, by acquiring multi-physical field parameters such as density field, velocity field, pressure field, temperature field, electron temperature field and magnetic field intensity, an instability sensitivity map containing instability intensity indicators and connection weights based on convective transmission, diffusion transmission and electromagnetic transmission mechanisms is constructed in the shock wave influence area, boundary layer area and their intersection area. Then, the map is subjected to multi-scale decomposition at the shock wave scale, vortex scale and molecular scale to extract the instability modes dominated by different physical mechanisms. The dominant instability modes are screened out using the hypersonic stability criterion considering the Mach number, Reynolds number and Hartmann number. Finally, an instability evolution prediction model including linear growth terms and nonlinear coupling terms is established to achieve accurate prediction of the spatiotemporal evolution law of hypersonic flow instability.

[0068] Specifically, in one embodiment of the present invention, step S1 includes:

[0069] First, obtain the density field ρ(x,y,z,t), which reflects the spatial variation of material distribution in the flow field. Second, obtain the velocity field The three components, streamwise, normal, and spanwise, characterize the flow's motion and shear properties. Third, the pressure field p(x, y, z, t) is obtained. This parameter is the basis for shock wave identification and pressure gradient calculation. Fourth, the temperature field T(x, y, z, t) is obtained, reflecting the thermodynamic state of the flow field.

[0070] In view of the characteristics of hypersonic multi-physics field coupling, it is also necessary to obtain the electron temperature field T e (x,y,z,t) and magnetic field strength The electron temperature field reflects the plasma effects produced by gas ionization under high temperature conditions, while the magnetic field intensity describes the interaction between the electromagnetic field and the flow field. All parameters are discretized using a unified spatial grid and time step to form a standardized multiphysics parameter dataset.

[0071] Specifically, in one embodiment of the present invention, step S2 includes:

[0072] S2.1 Basic structure of instability sensitivity map

[0073] Definition of graph: Instability sensitivity graph G = (V, E, W), where: V = v1, v2, ..., v n : A set of nodes, each node represents a spatial location in the flow field that is sensitive to instability; E = e ij : edge set, connecting node pairs with physical coupling relationship; W = w ij : A set of weights that quantifies the coupling strength between nodes.

[0074] Data structure representation: stored in the form of adjacency matrix, matrix A n×n Medium element A ij =w ij Represents the connection weight between node i and node j. If the two nodes are not connected, then A ij =0.

[0075] S2.2 Node settings and attribute definitions

[0076] By scanning the entire computational domain, nodes are set at locations that meet the following conditions:

[0077] Shock wave instability node: when When the shock wave instability node is set at the grid point (x, y, z), these nodes are mainly distributed in the pressure jump area before and after the shock wave. where ρ ∞ is the incoming flow density, U ∞ is the incoming flow velocity, L refThe threshold is determined based on the magnitude of the shock wave intensity to ensure that significant pressure jump areas can be captured.

[0078] Shear instability node: When When , set the shear instability nodes. These nodes are mainly distributed in the boundary layer and the shear layer area. The second preset threshold ε2 = 0.01 × U ∞ / δ, where δ is the boundary layer thickness. This threshold is determined based on the typical shear rate within the boundary layer, and is usually taken as 1% of the shear rate at the outer edge of the boundary layer as the identification standard.

[0079] Coupled instability node: It is set in the area where the above two conditions are met simultaneously (each takes a 50% threshold), representing the key position of the shock wave-boundary layer interaction.

[0080] Each node v i It has an instability strength index, and its calculation formula is:

[0081]

[0082] in, is the pressure gradient modulus, δ is the boundary layer thickness, p avg is the local average pressure, is the normal velocity gradient, U ∞ is the incoming flow velocity, T e is the electron temperature, T h is the heavy particle temperature.

[0083] The role of this indicator: Quantify the importance of the node: I node The larger the value, the stronger the instability sensitivity of the location, and the higher the weight in the subsequent modal analysis; modal extraction weight: during eigenvalue decomposition, the node instability intensity is used as the diagonal element of the weight matrix, affecting the calculation of the eigenvector; dominant mode screening: combined with the node instability intensity distribution, identify the mode that contributes most to the overall instability.

[0084] Two nodes v i and v j Conditions for establishing a connection:

[0085] Spatial distance criterion: |r i -r j | <R max , where R max =3δ is the maximum connection distance; Physical correlation criterion: at least one of the following conditions is met: Convective connection: there is significant convection transmission between two nodes (Pe ij >1); Diffusion connection: Diffusion transmission between two nodes cannot be ignored (|r i -r j |<2L diff); Electromagnetic connection: There is electromagnetic coupling effect (Re m,ij >0.1 and |r i -r j | <L em ).

[0086] Connection weight w ij Characterizes the physical coupling strength between two nodes:

[0087] w ij =W conv ·f conv (Pe ij ,|r i -r j |)+W diff ·f diff (Pe ij ,|r i -r j |)+W em

[0088] ·f em (Re m,ij ,|r i -r j |)

[0089] in, is the Peclet number, which characterizes the competition intensity between convection and diffusion, u i ,u j is the speed of nodes i and j, r i ,r j is the node position, v is the kinematic viscosity; is the magnetic Reynolds number, is the magnetic diffusivity, σ e is the electrical conductivity, μ0 is the vacuum magnetic permeability; W conv 、W diff 、W em are the weight coefficients of convective transmission, diffusion transmission and electromagnetic transmission respectively; f conv 、f diff 、f em are the corresponding convective transfer function, diffusion transfer function, and electromagnetic transfer function respectively. conv 、W diff 、W em The sum is 1.

[0090] The specific form of each transfer function is:

[0091]

[0092] is the convection characteristic length, Δt is the time step, and CFL is the dimensionless CFL number; is the diffusion characteristic length; is the characteristic length of plasma, c is the speed of light, ω p is the plasma frequency.

[0093] The algorithm flow for constructing the graph data structure is as follows:

[0094] Initialization: Create an empty node set V and adjacency matrix A; Node scanning: Traverse all grid points and determine the node position and type according to the threshold condition; Strength calculation: Calculate the instability strength index I for each node node ; Connection establishment: Check the connection criteria for all node pairs (i, j), and calculate the weight w if the conditions are met ij ; Matrix filling: Fill the weight values ​​into the adjacency matrix A to form the final graph data structure.

[0095] The constructed instability sensitivity map includes:

[0096] Node information table: Node coordinates (x i ,y i ,z i ), type (shock wave / shear / coupling), instability strength; adjacency matrix: n×n matrix, storing the connection weights between all nodes; degree matrix: diagonal matrix D, where D ii =∑ j w ij is the degree of node i.

[0097] Specifically, in one embodiment of the present invention, step S3 includes:

[0098] Hypersonic flow instabilities involve physical processes at multiple spatial and temporal scales. The shock wave scale primarily captures pressure wave propagation and shock wave structure changes, with its characteristic time scale being the acoustic time scale; the vortex scale primarily captures the evolution of rotating structures and shear layer instabilities, with its characteristic time scale being the convection time scale; and the molecular scale primarily captures viscous diffusion and heat conduction effects, with its characteristic time scale being the diffusion time scale. These three physical mechanisms are coupled under hypersonic conditions, collectively determining the complex behavior of flow instabilities.

[0099] To capture these physical processes at different scales, three targeted Laplacian operators are constructed. Each operator emphasizes the transmission characteristics of the corresponding scale by applying specific physical weighting to the basic connection weights of the instability sensitivity map.

[0100] The shock wave scale operator is mainly used to capture the propagation characteristics of pressure waves. Based on the basic connection weight matrix obtained in step S2, the weights are weighted according to the physical characteristics of the shock wave. The shock wave scale Laplace operator is:

[0101]

[0102] Shock wave weight calculation:

[0103]

[0104] where w ij is the basic connection weight of step S2, is the pressure gradient modulus at node i, is the reference pressure, L ref is the characteristic length of the object, |r i -r j | is the distance between nodes, L shock =L ref is the characteristic length of the shock wave. The pressure gradient intensity factor reflects the driving force strength of the shock wave propagation, and the distance attenuation factor reflects the spatial range of the shock wave influence. ref / L ref Determined based on the incoming flow pressure, the characteristic length of the shock wave is the characteristic size of the object, usually the chord length or diameter.

[0105] The vortex-scale operator is used to capture the evolution characteristics of the rotating structure and highlight the shear-driven instability mechanism.

[0106]

[0107] Vortex weight calculation:

[0108]

[0109] where ω i is the vorticity modulus length at node i, U ∞ is the incoming flow velocity, δ is the boundary layer thickness, L vortex =δ is the characteristic length of the vortex. The vorticity intensity factor is based on the vorticity at the node and the reference vorticity U ∞ The reference vorticity is determined by the ratio of δ / δ, which represents the typical shear rate at the outer edge of the boundary layer. The characteristic vortex length is the boundary layer thickness, usually defined as the 99% velocity boundary layer thickness, which reflects the main influence range of the vortex structure.

[0110] Molecular-scale operators are used to capture the viscous diffusion characteristics and reflect the effects of viscosity and heat conduction. Molecular-scale Laplace operator:

[0111]

[0112] Numerator weight calculation:

[0113]

[0114] Among them, Re δ,i =U i δ / v i is the local Reynolds number at node i, Ui is the velocity at the node, v i is the local kinematic viscosity, is the molecular characteristic length. The inverse of the Reynolds number reflects the relative importance of viscous effects, with viscous diffusion being more pronounced in the low Reynolds number region. The molecular characteristic length is determined based on viscous diffusion theory and reflects the characteristic influence range of viscous diffusion. Under hypersonic conditions, viscosity varies significantly due to temperature effects, necessitating consideration of the spatial distribution of local viscosity.

[0115] Perform eigenvalue decomposition on the three-scale Laplace operators:

[0116]

[0117] where scale∈{shock,vortex,molecular}, is the kth eigenvalue, is the corresponding eigenvector.

[0118] The physical meaning of eigenvalues ​​corresponds to the linear growth or decay rate of each mode. Positive eigenvalues ​​indicate unstable growth modes, while negative eigenvalues ​​indicate stable decay modes. The physical meaning of eigenvectors corresponds to the spatial distribution pattern of unstable perturbations, with each component of the vector representing the perturbation amplitude at the corresponding node.

[0119] In the numerical solution process, the Lanczos algorithm, which is specially designed for large sparse symmetric matrices, is used. This algorithm has good convergence and computational efficiency. The convergence criterion is set as the relative error of the eigenvalue is less than 10 -6 , the eigenvector residual is less than 10 -8 Mode number selection: Shock wave scale retains the first min (50, 0.1N total ) modes; the vortex scale retains the first min (100, 0.2N total ) modes; the molecular scale retains the first min(30,0.05N total ) modes. Among them, N total is the total number of graph nodes.

[0120] The modes at these three scales are combined to form a complete set of multiscale instability modes. Specifically, a complete modal matrix is ​​constructed, whose column vectors sequentially contain the modes at the shock wave scale, the vortex scale, and the molecular scale. Corresponding eigenvalue vectors and scale identifier vectors are also constructed. The former records the growth rate of each mode, and the latter identifies the physical scale to which each mode belongs (1-shock wave, 2-vortex, 3-molecule).

[0121] Specifically, in one embodiment of the present invention, step S4 includes:

[0122] Modal screening is performed using a stability criterion specifically tailored to hypersonic conditions. This criterion considers the combined effects of high Mach number, Reynolds number, and magnetic field effects on flow stability, and can accurately identify modes with the potential for instability growth.

[0123] Critical growth rate calculation:

[0124]

[0125] Among them, Ma is the Mach number, which represents the intensity of the compressibility effect; Re is the Reynolds number, which represents the ratio of inertial force to viscous force; Ha is the Hartmann number, which represents the intensity of the influence of the magnetic field on the flow.

[0126] Parameter calculation method: Mach number Ma=U ∞ / a ∞ , where U ∞ is the incoming flow velocity, a ∞ is the incoming sound velocity; Reynolds number Re=ρ ∞ U ∞ L ref / μ ∞ , where ρ ∞ is the incoming flow density, L ref is the reference length, μ ∞ is the dynamic viscosity of the incoming flow; Hartmann number Where B is the magnetic induction intensity, σ e is the conductivity.

[0127] The critical growth rate is the threshold at which the mode begins to grow unstably under given hypersonic flow conditions. Only when the eigenvalue of the mode is greater than this critical value does the mode have the potential to grow unstably.

[0128] All modes (including modes of three scales) are judged for their critical growth rate and stable modes that do not meet the instability conditions are eliminated. Screening condition: λ k >λ crit ; where λ k is the eigenvalue of the kth mode, i.e., the linear growth rate of the mode. Under hypersonic conditions, many modes remain stable due to strong viscous dissipation or magnetic field stabilization, and need to be identified and eliminated by this criterion.

[0129] After judging by the critical growth rate, the candidate unstable mode set φ is obtained k , in is the set of modal indices that meet the instability conditions.

[0130] For the candidate unstable modes determined by the critical growth rate, a comprehensive importance evaluation is performed based on three dimensions: modal growth rate, spatial influence range, and energy contribution.

[0131] The modal growth rate directly reflects the growth rate of the unstable mode and is the primary indicator for judging the importance of the mode.

[0132]

[0133] Among them S growth,k The growth rate score of the kth mode is in the range of (0,1). The higher the growth rate of the mode, the closer its score is to 1, indicating that the mode has a stronger ability to grow in an unstable manner in the linear stage.

[0134] The spatial influence range measures the distribution breadth of the instability mode in the physical space. The mode with a large influence range has a more significant disturbing effect on the overall flow field. Calculation of spatial influence range:

[0135]

[0136] in is the amplitude of the kth mode at node i, r i is the distance from node i to the geometric center, N total is the total number of nodes. Normalized score of influence range:

[0137]

[0138] This index reflects the spatial distribution characteristics of the mode. The mode with high localization degree R k Smaller value and higher degree of globalization of modal R k The value is larger.

[0139] Energy contribution evaluation modal energy contribution in the entire instability process, weighted calculation based on the node instability intensity. Modal energy contribution calculation:

[0140]

[0141] in is the node instability strength weight matrix, I i is the instability strength index of node i. Energy contribution E k It directly reflects the relative importance of the mode in the instability process. Its value range is (0,1), and the sum of the energy contributions of all modes is 1.

[0142] The three evaluation dimensions are weighted and combined to construct a comprehensive importance index. Comprehensive importance calculation:

[0143] S k =α·S growth,k +β·S range,k +γ·E k

[0144] Where α, β, and γ are weight coefficients, satisfying α+β+γ=1.

[0145] Sort the candidate unstable modes based on the comprehensive importance index and select the first several modes as the dominant unstable modes. Mode sorting: Sort all candidate modes according to the comprehensive importance S k Sort from large to small to get an ordered mode sequence. Determine the number of dominant modes: Threshold method: Select S k >S threshold All modes of S threshold =0.1, which is suitable for the case where the modal importance is relatively evenly distributed; cumulative contribution method: select the modes whose cumulative comprehensive importance reaches 80% of the total to ensure that the dominant mode set covers the main instability characteristics.

[0146] The output of step S4 includes the set of dominant instability modes and their related information. Dominant mode matrix: Each column is the spatial distribution of a dominant instability mode. Dominant mode eigenvalue: Corresponding to the linear growth rate of each dominant mode. Comprehensive importance vector: Record the comprehensive importance score of each dominant mode. Scale identification vector: Identify the physical scale to which each dominant mode belongs (1-shock wave, 2-vortex, 3-molecule).

[0147] Specifically, in one embodiment of the present invention, step S5 includes:

[0148] The evolution process of hypersonic flow instability can be described by the modal decomposition method. In the linear stage, each unstable mode evolves independently, and its growth rate is determined by the corresponding eigenvalue. As the instability develops, nonlinear interactions begin to occur between the modes, which determine the late behavior of the instability evolution and the final saturation state. Based on the modal decomposition theory, any instability field variable can be expressed as a linear superposition of the dominant instability modes. The time evolution of each mode is controlled by an ordinary differential equation, which contains a linear growth term and a nonlinear coupling term. By solving this set of coupled evolution equations, the evolution law of the instability field variables over time can be obtained. The instability evolution prediction model adopts the form of an evolution equation based on the linear superposition of the dominant modes:

[0149]

[0150] Among them, k is the time-dependent amplitude of the kth dominant instability mode, indicating the intensity change of the mode during the evolution process; k is the linear growth rate, i.e. the dominant modal eigenvalue selected in step S4; N kmn is the nonlinear coupling coefficient, which describes the nonlinear influence of the mth and nth modes on the evolution of the kth mode.

[0151] In this equation, the evolution of each mode contains two contributions. The first part is λ k Ψ k is a linear term that describes the self-growth or attenuation of the mode; the second part ∑ m,n N kmn Ψ m Ψ n is a nonlinear term that describes the effect of the interaction between different modes on the evolution of the mode.

[0152] Nonlinear coupling coefficient N kmn The calculation is based on the physical mechanism of hypersonic flow and the spatial distribution characteristics of the modes:

[0153]

[0154] where φ k (r) is the spatial distribution function of the kth dominant mode, is a nonlinear operator, Ω is the computational domain, and r is the spatial position vector.

[0155] Nonlinear operators The nonlinear terms based on the Navier-Stokes equations are constructed as follows:

[0156]

[0157] where u m ,u n are the velocity disturbances corresponding to the mth and nth modes, respectively, p m ,p n is the corresponding pressure disturbance, and ρ0 is the reference density.

[0158] In order to solve the instability evolution equation, it is necessary to set appropriate initial conditions and boundary conditions.

[0159] Initial condition setting: in is the initial amplitude of the kth mode. The determination of the initial amplitude is based on the initial disturbance distribution in the flow field: in is the initial perturbation velocity field distribution.

[0160] Methods for setting initial disturbances: Random disturbance method: adding small random disturbances to the flow field, the amplitude of which is usually 0.1%-1% of the incoming flow velocity; specific disturbance method: setting a specific form of initial disturbance based on actual physical problems, such as disturbances caused by wall roughness; experimental data method: setting initial conditions based on disturbance data measured experimentally.

[0161] The fourth-order Runge-Kutta method is used to solve the evolution equations.

[0162] Time step determination:

[0163]

[0164] Among them C CFL =0.1-0.5 is the CFL coefficient, λ max is the maximum linear growth rate, N max is the maximum nonlinear coupling coefficient, |Ψ max | is the maximum modal amplitude at the current moment.

[0165] Convergence criterion: When the relative change rate of all modal amplitudes is less than the set threshold, the calculation is considered to have converged:

[0166]

[0167] where ∈ conv =10 -6 is the convergence threshold, and the superscript n represents the number of time steps.

[0168] After solving the evolution equations to obtain the time-dependent amplitudes of each dominant mode, the complete spatiotemporal distribution of the instability field variables can be reconstructed.

[0169] Field variable reconstruction:

[0170]

[0171] where u disturb (r, t) is the spatial and temporal distribution of the disturbance velocity field, p disturb (r, t) is the spatial and temporal distribution of the disturbance pressure field, p k (r) is the pressure distribution corresponding to the kth mode, N dominant is the total number of dominant modes.

[0172] Full field variable calculation:

[0173] u(r,t)=u base (r)+u disturb (r,t)

[0174] p(r,t)=p base (r)+p disturb (r,t)

[0175] where u base (r) and p base (r) are the velocity and pressure distributions of the basic flow field, respectively.

[0176] The output of step S5 includes the spatiotemporal distribution prediction results of instability evolution and related analysis indicators. Time series output: Output the evolution curve Ψ of each dominant mode amplitude over timek (t), used to analyze the growth and attenuation trends of each mode. Spatial distribution output: output the spatial distribution u of the disturbance field at different times disturb (r,t i ) and p disturb (r,t i ), where t i is the selected time series. Identification of the unstable development stage: Linear stage: When |∑ m,n N kmn Ψ m Ψ n |<0.1|λ k Ψ k |, the mode is in the linear development stage; nonlinear stage: when the nonlinear term begins to dominate, the instability enters the nonlinear development stage; saturation stage: when When , the mode reaches saturation.

[0177] like Figure 3 As shown, the present invention also provides a hypersonic flow instability multi-field coupling prediction system, which is used to implement any of the above methods, and the system includes:

[0178] A data acquisition module is used to acquire and pre-process the multi-physics parameters of the hypersonic flow field;

[0179] A graph construction module is used to set graph nodes in the shock wave impact area, boundary layer area, and their intersection area according to preset thresholds, calculate the node instability intensity index, and calculate the connection weights between nodes based on convective transmission, diffusion transmission, and electromagnetic transmission mechanisms to construct the instability sensitivity graph data structure;

[0180] The multi-scale decomposition module is used to construct feature extraction operators at the shock wave scale, vortex scale, and molecular scale, perform eigenvalue decomposition on each scale operator, extract the corresponding unstable modes, and merge them to form a complete unstable mode set;

[0181] Mode screening module, used to determine the dominant instability mode based on critical growth rate determination and comprehensive importance evaluation;

[0182] The prediction module is used to construct an instability evolution prediction model including nonlinear interactions based on the dominant instability mode, and solve the evolution equation to output the spatiotemporal distribution law of instability.

[0183] It should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit the present invention. Although the present invention has been described in detail with reference to the preferred embodiments, those skilled in the art should understand that the technical solutions of the present invention may be modified or replaced by equivalents without departing from the spirit and scope of the technical solutions of the present invention, which should all be included in the scope of the claims of the present invention.

Claims

1. A multi-field coupling prediction method for hypersonic flow instability, characterized by: include: S1. Obtaining the multi-physics parameter distribution of the hypersonic flow field; S2. Based on the shock wave-boundary layer interaction mechanism, the instability sensitivity map data structure is constructed with multi-physics field parameters as input. Graph nodes are set in the shock wave impact area, boundary layer area, and their intersection area, and the connection weights between nodes are calculated. S3. Perform multi-scale decomposition on the instability sensitivity map data structure to extract the instability modes at the shock wave scale, vortex scale, and molecular scale; S4. Use hypersonic stability criteria to screen the instability modes and determine the dominant instability mode; S5. Construct an instability evolution prediction model based on the dominant instability mode and output the spatiotemporal evolution law of hypersonic flow instability.

2. The hypersonic flow instability multi-field coupling prediction method according to claim 1, characterized in that: The method for setting the graph node in step S2 includes: Setting a shock wave instability node in an area where the pressure gradient exceeds a first preset threshold; Setting a shear instability node in an area where the velocity gradient exceeds a second preset threshold; Set up a coupling instability node in the intersection area of ​​shock wave and boundary layer; The instability strength index of each node is calculated according to the following formula: in, is the pressure gradient modulus, δ is the boundary layer thickness, p avg is the local average pressure, is the normal velocity gradient, U ∞ is the incoming flow velocity, T e is the electron temperature, T h is the heavy particle temperature.

3. The multi-field coupling prediction method for hypersonic flow instability according to claim 1, characterized in that: The connection weights between nodes in step S2 are determined based on the physical transmission mechanism, and the weight calculation formula is: w ij =W conv f conv (Per ij , i -r j |)+W diff f diff (Per ij , i -r j |)+W em ·f em (Re m,ij ,|r i -r j |) in, is the Peclet number, which characterizes the competition intensity between convection and diffusion, u i ,u j is the speed of nodes i and j, r i ,r j is the node position, v is the kinematic viscosity; is the magnetic Reynolds number, is the magnetic diffusivity, σ e is the electrical conductivity, μ0 is the vacuum magnetic permeability; W conv 、W diff 、W em are the weight coefficients of convective transmission, diffusion transmission and electromagnetic transmission respectively; f conv 、f diff 、f em are the corresponding convective transfer function, diffusion transfer function, and electromagnetic transfer function, respectively.

4. The multi-field coupling prediction method for hypersonic flow instability according to claim 3, characterized in that: The specific form of each transfer function is: is the convection characteristic length, Δt is the time step, and CFL is the dimensionless CFL number; is the diffusion characteristic length; is the characteristic length of plasma, c is the speed of light, ω p is the plasma frequency.

5. The hypersonic flow instability multi-field coupling prediction method according to claim 1, characterized in that: The multi-scale decomposition in step S3 includes: Construct Laplace operators at three scale levels: shock-wave scale operators are used to capture pressure wave propagation characteristics, vortex-scale operators are used to capture the evolution characteristics of rotating structures, and molecular-scale operators are used to capture viscous diffusion characteristics. Perform eigenvalue decomposition on the Laplace operator of each scale to obtain the unstable mode of the corresponding scale; The modes at the three scales are merged to form a complete multi-scale unstable mode set.

6. The multi-field coupling prediction method for hypersonic flow instability according to claim 1, characterized in that: The hypersonic stability criterion in step S4 includes a critical growth rate determination, and the critical growth rate is calculated as follows: Among them, Ma is the Mach number, Re is the Reynolds number, and Ha is the Hartmann number.

7. The multi-field coupling prediction method for hypersonic flow instability according to claim 1, characterized in that: The determination of the dominant instability mode in step S4 is based on a comprehensive evaluation of the modal growth rate, spatial influence range and energy contribution.

8. The multi-field coupling prediction method for hypersonic flow instability according to claim 1, characterized in that: The instability evolution prediction model in step S5 includes: The instability evolution equation based on the linear superposition of the dominant modes is established in the following form: Among them, ψ k is the time-dependent amplitude of the kth dominant instability mode, λ k is the linear growth rate, N kmn is the nonlinear coupling coefficient; The spatiotemporal distribution prediction results of the instability field variables are obtained by solving the evolution equation.

9. The multi-field coupling prediction method for hypersonic flow instability according to claim 1, characterized in that: The multi-physics field parameters in step S1 include density field, velocity field, pressure field, temperature field, electron temperature field and magnetic field intensity.

10. Hypersonic flow instability multi-field coupling prediction system, characterized by: The system is used to implement the method according to any one of claims 1 to 9, and the system includes: A data acquisition module is used to acquire and pre-process the multi-physics parameters of the hypersonic flow field; A graph construction module is used to set graph nodes in the shock wave impact area, boundary layer area, and their intersection area according to preset thresholds, calculate the node instability intensity index, and calculate the connection weights between nodes based on convective transmission, diffusion transmission, and electromagnetic transmission mechanisms to construct the instability sensitivity graph data structure; The multi-scale decomposition module is used to construct feature extraction operators at the shock wave scale, vortex scale, and molecular scale, perform eigenvalue decomposition on each scale operator, extract the corresponding unstable modes, and merge them to form a complete unstable mode set; Mode screening module, used to determine the dominant instability mode based on critical growth rate determination and comprehensive importance evaluation; The prediction module is used to construct an instability evolution prediction model including nonlinear interactions based on the dominant instability mode, and solve the evolution equation to output the spatiotemporal distribution law of instability.

Images