Modeling method of shock wave light field in hypersonic flow field

Abstract

The application belongs to the field of numerical simulation, and particularly relates to a modeling method of a shock wave light field in a hypersonic flow field, comprising the following steps: adopting an unsupervised learning algorithm to perform gradient clustering on flow field data, and determining the position of a shock wave and the inner and outer boundaries of a shock wave domain in the flow field domain; constructing a hybrid dimension numerical continuous and discrete strategy, determining the main direction of a light ray based on the strategy, performing continuous medium approximation processing on the main direction of the light ray, and performing structured grid discrete processing on the non-main direction of the light ray; based on the analogy of Fermat's principle and Maupertuis' principle, the transmission of the light ray in the refractive index field is analogized to the movement of equivalent real particles in the refractive index potential field, and the dynamic relationship between the equivalent real particles and the refractive index potential field is established; based on the dynamic relationship, the transmission path of the light ray in the shock wave is simulated, and the modeling of the shock wave light field in the hypersonic flow field is realized. The application significantly improves the efficiency and accuracy of the modeling of the shock wave light field in the hypersonic flow field.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of numerical simulation, and particularly relates to a modeling method of a shock wave light field in a hypersonic flow field. BACKGROUND

[0002] With the rapid development of aircraft technology, the flight speed has been increased from subsonic to hypersonic range. In the hypersonic flight state, the flow field around the aircraft presents strong space-time dynamic characteristics, resulting in deflection and high-frequency jitter of the light beam passing through the flow field, and then causing imaging blur and jitter, i.e. aerodynamic optical effect.

[0003] In order to study the aerodynamic optical effect, the traditional wind tunnel experiment method is adopted, but this method is high in cost and difficult to flexibly adjust the space-time parameters. Therefore, the numerical simulation method becomes an important means for studying the light field transmission behavior in the hypersonic flow field.

[0004] In the prior art, the algorithm based on light ray tracing is widely used in light field modeling. This kind of method solves the light ray equation by high-order Runge-Kutta method, and needs to repeatedly interpolate the flow field data and its density gradient, which is large in calculation resource consumption and low in efficiency. Especially in the hypersonic flow field, in order to ensure the space-time resolution, the computational fluid dynamics simulation needs to adopt nanosecond time step, resulting in huge data amount, and the traditional light ray tracing method is difficult to meet the real-time and high efficiency requirements of dynamic light field modeling.

[0005] In addition, although high-precision flow field simulation can generate ten million grid data, not all flow field regions have significant influence on optical transmission. Among them, the shock wave region plays a leading role in light deflection. If the light ray tracing is performed point by point in the whole flow field, a large amount of redundant calculation will be caused, and the rapid modeling of the light field cannot be realized.

[0006] Therefore, there is an urgent need for a method capable of efficiently and accurately modeling the light field in the shock wave region in the hypersonic flow field to overcome the shortcomings of the prior art. SUMMARY

[0007] Therefore, the present application aims to provide a modeling method of a shock wave light field in a hypersonic flow field, which can efficiently and accurately model the light field in the shock wave region in the hypersonic flow field to overcome the shortcomings of the prior art.

[0008] To achieve the above-mentioned purpose, the technical scheme of the present application is as follows:

[0009] A modeling method of a shock wave light field in a hypersonic flow field, comprising the following steps:

[0010] S1: using an unsupervised learning algorithm to perform gradient clustering on the flow field data to determine the position of the shock wave in the flow field domain, the inner and outer boundaries of the shock wave domain, and to obtain the inner and outer boundaries of the shock wave domain and the spatial distribution analytical expression of the shock wave by numerical fitting;

[0011] S2: constructing a hybrid dimension numerical continuous and discrete strategy, determining the main direction of the light ray based on the strategy, performing continuous medium approximation processing on the main direction of the light ray, and performing structured grid discrete processing on the non-main direction of the light ray;

[0012] S3: based on the analogy of Fermat's principle and the Mopaidu principle under the condition of energy conservation, the transmission of light in the refractive index field is analogous to the motion of equivalent real particles in the refractive index potential field, thereby establishing the dynamic relationship between the equivalent real particles and the refractive index potential field;

[0013] S4: based on the dynamic relationship between the equivalent real particles and the refractive index potential field, the spatial distribution analytical expression of the inner boundary of the shock wave domain is solved by combining the straight line equation of the light ray and the shock wave domain, the transmission path of the light ray in the shock wave is simulated, and the modeling of the shock wave light field in the hypersonic flow field is realized.

[0014] Further, in step S1, the specific process of using an unsupervised learning algorithm to perform gradient clustering on the flow field data is as follows:

[0015] S11: abstracting and simplifying the flow field data into a one-dimensional data set; wherein the flow field data is composed of node coordinates and flow field attribute data, and the node position is constant;

[0016] S12: randomly selecting two nodes in the flow field data as initial centroids, calculating the density value of each node in the flow field data and the Euclidean distance between the two initial centroids, and according to the size of the Euclidean distance, classifying the density value of each node into the cluster corresponding to the nearest initial centroid;

[0017] S13: recalculating the new centroid of each cluster to replace the old centroid, repeating steps S11 and S12 until the centroid and cluster assignment converge and no longer change, and dividing the flow field domain into far field and near field;

[0018] S14: taking the demarcation between the far field and the near field as the outer boundary of the shock wave domain, performing a downward translation operation on the outer boundary of the shock wave domain, and when the refractive index gradient reaches a peak value, taking the position of the peak value as the position of the shock wave;

[0019] S15: selecting a node in the outer boundary of the shock wave domain symmetric to the position of the shock wave as the inner boundary of the shock wave domain.

[0020] Further, in step S2, the determination method of the main direction of the light ray is as follows:

[0021] In the two-dimensional flow field, the angles between the light ray and the x-axis, y-axis are compared, and the coordinate axis with the smallest angle is selected as the main direction of the light ray;

[0022] In the three-dimensional flow field, the angles between the light ray and the x-axis, y-axis, z-axis are compared, and the coordinate axis with the smallest angle is selected as the main direction of the light ray.

[0023] Further, in step S2, the specific process of implementing the continuous medium approximation processing on the main direction of the light ray is as follows:

[0024] A point is selected on the outer boundary, shock wave, and inner boundary of the shock wave domain along the main direction of the light ray, and a polynomial equation about the main direction coordinate is constructed according to the coordinates and refractive index values of the three points. The polynomial equation is solved by the Cramer rule to obtain an analytical expression of the continuous distribution of the refractive index of the shock wave domain in the main direction of the light ray, thereby realizing the continuous medium approximation processing on the main direction of the light ray.

[0025] Further, in step S2, the specific process of implementing the structured grid discrete processing on the non-main direction of the light ray is as follows:

[0026] In the two-dimensional flow field, the non-main direction of the light ray is divided into equal-interval micro-intervals, and the flow field in a single micro-interval remains continuous in the main direction of the light ray.

[0027] In the three-dimensional flow field, the non-main direction of the light ray is divided into micro-units of equal size, and the flow field in a single micro-unit remains continuous in the main direction of the light ray.

[0028] Further, in step S3, the specific process of establishing the dynamic relationship between the equivalent real particles and the refractive index potential field is as follows:

[0029] The variational form of Fermat's principle is:

[0030] ;

[0031] wherein, is the refractive index, is the light path;

[0032] The Maupertuis principle under the energy conservation condition is:

[0033] ;

[0034] wherein, is the optical momentum, is the generalized coordinate;

[0035] By analogy between Fermat's principle and Maupertuis' principle under the condition of energy conservation, the transmission of light in the refractive index field is analogous to the motion of an equivalent real particle in the refractive index potential field, and the mass of the equivalent real particle is set to 1, so the dynamic relationship between the equivalent real particle and the refractive index potential field is:

[0036] ;

[0037] wherein, is the velocity of the equivalent real particle.

[0038] Further, in step S4, the refractive index of the initial intersection point is set to , and the specific process of simulating the transmission path of the light ray in the shock domain is:

[0039] If the flow field is a two-dimensional flow field and the main direction of the light ray is the y-axis, the velocity component expressions of the equivalent real particle in the main direction and the non-main direction of the light ray are obtained based on the dynamic relationship between the equivalent real particle and the refractive index potential field:

[0040] ;

[0041] ;

[0042] wherein, is the velocity component of the equivalent real particle in the non-main direction of the light ray, is the velocity component of the equivalent real particle in the main direction of the light ray, is the refractive index distribution along the main direction of the light ray in the two-dimensional flow field, is the angle between the light ray at the initial intersection point and the x-axis;

[0043] The position differential expression of the equivalent real particle is derived from the velocity component expressions of the equivalent real particle in the main direction and the non-main direction of the light ray:

[0044] ;

[0045] The position expression is obtained by integrating the position differential expression of the equivalent real particle on both sides:

[0046] ;

[0047] When the equivalent real particle passes through a small interval in the non-main direction of the light ray, the refractive index n in the small interval is updated, and the motion trajectory of the equivalent real particle in the small interval is calculated based on the updated refractive index n;

[0048] If the flow field is a three-dimensional flow field and the main direction of the light ray is the y-axis, the velocity component expressions of the equivalent physical particle in the non-main direction plane and the main direction of the light ray are obtained based on the dynamic relationship between the equivalent physical particle and the refractive index potential field:

[0049] ;

[0050] ;

[0051] wherein, is the velocity component of the equivalent physical particle in the x-z plane, is the velocity component of the equivalent physical particle in the main direction of the light ray, is the refractive index distribution along the main direction of the light ray in the three-dimensional flow field, is the angle between the light ray and the x-z plane at the initial intersection point;

[0052] The position differential expression of the equivalent physical particle is derived based on the velocity component expressions of the equivalent physical particle in the x-z plane and the y-axis:

[0053] ;

[0054] wherein, is the projection of the equivalent physical particle in the x-z plane;

[0055] The position expression is obtained by integrating the position differential expression of the equivalent physical particle:

[0056] ;

[0057] When the equivalent physical particle passes through a micro-unit in the non-main direction of the light ray, the refractive index n in the micro-unit is updated, and the motion trajectory of the equivalent physical particle in the micro-unit is calculated based on the updated refractive index n.

[0058] Further, in step S4, the specific process of solving the initial intersection point of the light ray and the shock domain is that the straight line equation of the light ray and the spatial distribution analytical expression of the inner boundary of the shock domain are solved together:

[0059] In a two-dimensional flow field, the straight line equation of the light ray is converted into , and the spatial distribution analytical expression of the inner boundary of the shock domain is converted into , and the initial intersection point of the straight line propagation of the light ray and the inner boundary of the shock domain is obtained by solving the two together:

[0060] ;

[0061] In a three-dimensional flow field, the straight line equation of the light ray is converted into , and the spatial distribution analytical expression of the inner boundary of the shock domain is converted into , both together can get the initial intersection point of the light along the straight line propagation and the inner boundary of the shock domain:

[0062] .

[0063] Further, in step S4, the specific process of modeling the shock light field in the hypersonic flow field is:

[0064] Calculate the optical path length of each light , the calculation formula is:

[0065] ;

[0066] Wherein, s i is the propagation distance of the light in each tracking step, n i is the refractive index of each tracking step;

[0067] Based on the optical path length of each light, calculate the optical path difference of each light , the calculation formula is:

[0068] ;

[0069] Wherein, is the average value of all ;

[0070] Based on the optical path difference of each light , calculate the wavefront aberration , the calculation formula is:

[0071] ;

[0072] Wherein, k = 2π / λ , λ is the wavelength of the light;

[0073] Based on the wavefront aberration , calculate the actual pupil function , the calculation formula is:

[0074] ;

[0075] Wherein, is the ideal pupil function, indicates the phase disturbance factor;

[0076] Based on the actual pupil function , calculate the complex amplitude function , the calculation formula is:

[0077] ;

[0078] Based on the complex amplitude function , the point spread function is calculated , the calculation formula is:

[0079] ;

[0080] Wherein, is the complex conjugate of ;

[0081] Based on the point spread function , the optical transfer function is calculated , the calculation formula is:

[0082] ;

[0083] Based on the optical transfer function , the modulation transfer function is calculated , the calculation formula is:

[0084] .

[0085] Further, the transmission path of multiple light rays in the shock wave is simulated, and the modeling of the shock wave light field in the hypersonic flow field is completed based on the transmission path of multiple light rays.

[0086] Compared with the prior art, the application can achieve the following beneficial effects:

[0087] (1) The calculation efficiency is significantly improved: the present application focuses on the shock wave area which plays a leading role in light field transmission through shock wave domain identification and extraction technology, avoiding the problem of redundant calculation of the entire large flow field in traditional methods. Combined with the mixed dimension numerical strategy (continuous light ray main direction, discrete light ray non-main direction), the number of complex interpolation operations is greatly reduced, thereby fundamentally reducing the calculation complexity, realizing the rapid modeling of the light field of the hypersonic flow field, and being especially suitable for dynamic process analysis.

[0088] (2) High modeling accuracy: based on the gradient clustering method of unsupervised learning, the shock wave interface and range can be accurately identified in a data-driven manner, overcoming the limitations of manually setting thresholds. At the same time, the Hamiltonian optical principle is introduced, and the light wave propagation behavior in the time-varying refractive index medium is equivalent to the dynamic process of particles with mass in the potential field, providing a rigorous dynamic theoretical basis for light field tracing, and improving the modeling accuracy of the shock wave light field.

[0089] (3) Low resource consumption, strong practicability: the present application greatly reduces the demand for computer memory and computing resources under the premise of guaranteeing spatial resolution by focusing on the key area (shock wave domain) and adopting an efficient mixed dimension strategy, which makes it possible to realize high-resolution flow field light field modeling of the order of ten million grids on ordinary computing platforms, reduces the technical threshold and application cost, and has high engineering application value.

[0090] (4) Provide an effective tool for dynamic light field analysis: the traditional method is difficult to process high time resolution dynamic flow field sequence due to low computational efficiency. The efficient characteristics of the present application enable it to seamlessly connect the nanosecond time step data generated by computational fluid dynamics simulation, thereby realizing effective simulation and analysis of the dynamic, high-frequency evolution process of the hypersonic flow field aerodynamic optical effect (such as image blur, jitter), and providing a powerful numerical simulation tool for aircraft optical guidance, observation system performance evaluation, etc. BRIEF DESCRIPTION OF DRAWINGS

[0091] The accompanying drawings, which form a part of the present application, are used to provide further understanding of the present application, and the illustrative embodiments of the present application and their description are used to explain the present application, and do not constitute improper limitations on the present application. In the drawings:

[0092] Figure 1 The flowchart of the modeling method of the shock wave light field in the hypersonic flow field according to the embodiment of the present application is shown in the figure;

[0093] Figure 2 The refractive index distribution diagram of the flow field domain according to the embodiment of the present application is shown in the figure;

[0094] Figure 3 The refractive index gradient distribution diagram of the flow field domain according to the embodiment of the present application is shown in the figure;

[0095] Figure 4 The demarcation diagram of the near field and the far field according to the embodiment of the present application is shown in the figure;

[0096] Figure 5 The shock wave domain diagram according to the embodiment of the present application is shown in the figure;

[0097] Figure 6 The transmission path diagram of the simulated light in the shock wave in the two-dimensional flow field according to the embodiment of the present application is shown in the figure;

[0098] Figure 7 The transmission path diagram of the simulated light in the shock wave in the three-dimensional flow field according to the embodiment of the present application is shown in the figure. DETAILED DESCRIPTION

[0099] In order to make the purposes, technical solutions and advantages of the present application clearer, the present application will be further described in detail below in combination with the drawings and specific embodiments. It should be understood that the specific embodiments described herein are only used to explain the present application, and do not constitute a limitation on the present application.

[0100] It should be noted that the embodiments in the present application and the features in the embodiments can be combined with each other without conflict.

[0101] In the description of the present application, it should be understood that the terms "center", "longitudinal", "transverse", "upper", "lower", "front", "rear", "left", "right", "vertical", "horizontal", "top", "bottom", "inner", "outer" and the like indicate the orientation or positional relationship shown in the drawings, and are only for the convenience of describing the present application and simplifying the description, and do not indicate or imply that the device or element referred to must have a particular orientation, be constructed and operated in a particular orientation, and therefore cannot be understood as a limitation on the present application. In addition, the terms "first", "second" and the like are only for descriptive purposes and cannot be understood as indicating or implying relative importance or implicitly indicating the number of technical features indicated. Therefore, the features defined with "first", "second" and the like can explicitly or implicitly include one or more of the features. In the description of the present application, unless otherwise specified, the meaning of "a plurality of" is two or more.

[0102] In the description of the present application, it should be noted that unless otherwise specified and limited, the terms "assembly", "connection", "connection" should be understood broadly, for example, it can be fixedly connected, or it can be detachably connected, or integrally connected; it can be mechanically connected, or it can be electrically connected; it can be directly connected, or it can be indirectly connected through an intermediate medium, or it can be connected inside two elements. For those skilled in the art, the specific meaning of the above terms in the present application can be understood according to the specific circumstances.

[0103] The present application will be described in detail below with reference to the drawings and in combination with the embodiments.

[0104] As Figure 1 shown, the present application provides a modeling method of shock wave light field in hypersonic flow field, comprising the following steps:

[0105] S1: using an unsupervised learning algorithm to perform gradient clustering on the flow field data, determining the position of the shock wave in the flow field domain, the inner and outer boundaries of the shock wave domain, and obtaining the inner and outer boundaries of the shock wave domain and the spatial distribution analytical expression of the shock wave by numerical fitting.

[0106] The present application adopts an unsupervised learning algorithm to perform gradient clustering on flow field data, establishes a shock interface identification criterion taking the maximum density gradient modulus as the criterion, performs clustering analysis on the flow field data according to the refractive index gradient of the flow field, and further accurately identifies and separates the shock wave with the most significant influence on the light field in the flow field domain, determines the inner boundary, outer boundary and position of the shock wave domain, and obtains the spatial distribution analytical expression of the inner boundary, outer boundary and shock wave of the shock wave domain through numerical fitting. The refractive index distribution and refractive index gradient distribution of the flow field domain are as shown in Figs. 2 and Figure 3 .

[0107] The specific process of adopting an unsupervised learning algorithm to perform gradient clustering on flow field data is as follows:

[0108] S11: Abstract and simplify the flow field data into a one-dimensional data set.

[0109] The flow field data is composed of node coordinates (used to describe the position information of the flow field nodes in two-dimensional / three-dimensional space) and flow field attribute data (used to reflect the core parameters of the physical state of the flow field, such as density, refractive index, etc.), which belongs to a multi-dimensional data set. Under the premise that the node position remains constant, the core attribute values (for example, density values) of each node are extracted in a fixed order (such as the node sequence along a certain coordinate axis direction or the node number order of the flow field data), and these core attribute values are sequentially arranged according to a predetermined rule. The multi-dimensional flow field data is abstracted and simplified into a one-dimensional data set that only reflects the continuous change of the core attribute values, which can be directly used for subsequent Euclidean distance calculation, greatly reducing the clustering complexity. The flow field domain covers fluids in different flow states, and the refractive index of fluids in different states differs significantly. Based on this, the fluids in the flow field domain can be roughly divided into far field and near field.

[0110] S12: Randomly select two nodes in the flow field data as initial centroids, calculate the Euclidean distance between the density value of each node in the flow field data and the two initial centroids, and according to the size of the Euclidean distance, classify the density value of each node into the cluster corresponding to the initial centroid closest to it.

[0111] In a one-dimensional space, the Euclidean distance d of the density value of a certain node and the initial centroid is calculated according to the following formula:

[0112] .

[0113] S13: Recalculate the new centroid of each cluster to replace the old centroid, and repeat the operations of density value reassignment and centroid update until the assignment of the centroid and the cluster converges and no longer changes, and the flow field domain is divided into far field and near field, and the demarcation between the near field and the far field is determined.

[0114] The demarcation between the near field and the far field is as shown in​​​Figure 4 As shown.

[0115] S14: Taking the demarcation between the far field and the near field as the outer boundary of the shock domain, performing a downward translation operation on the outer boundary of the shock domain, and taking the position of the peak value when the refractive index gradient reaches the peak value as the position of the shock wave.

[0116] S15: Selecting a node symmetric to the position of the shock wave in the outer boundary of the shock domain as the inner boundary of the shock domain.

[0117] At this point, the shock domain is completely separated, as shown in the figure, thereby achieving effective identification and extraction of the shock wave. The shock wave has a certain form in the hypersonic flow field. Figure 5 The inner boundary of the shock domain, the outer boundary of the shock domain and the spatial distribution of the shock wave can be obtained by numerical fitting. Specifically, the shock wave position point set (the point with the maximum refractive index gradient), the outer boundary point set of the shock domain and the inner boundary point set of the shock domain are extracted from the flow field data, and a suitable analytical function form is selected for fitting according to the physical characteristics (such as symmetry, curvature, etc.) of the shock wave form, for example: polynomial function (such as quadratic, cubic), exponential function, Gaussian function, spline function.

[0118] The present application focuses on the shock wave region which plays a leading role in the light field transmission by identifying and extracting the shock wave domain, avoiding the problem of redundant calculation of the entire large flow field in the traditional method.

[0119] S2: Construct a mixed dimension numerical continuous and discrete strategy, determine the main direction of the light based on the strategy, implement continuous medium approximation processing on the main direction of the light, and implement structured grid discrete processing on the non-main direction of the light.

[0120] In optical flow field analysis, when using the reverse optical tracing method, the light starts from the entrance pupil of the optical system and is projected into the flow field.

[0121] In a two-dimensional flow field, the exit angle of the light is

[0122] . Wherein, is the angle between the light and the x-axis, is the angle between the light and the y-axis. If , the y-axis direction is the main direction of the light; if , the x-axis direction is the main direction of the light.

[0123] In a three-dimensional flow field, the exit angle of the light is . Wherein, is the angle between the light and the x-axis, is the angle between the light and the y-axis, ​Let θ be the angle between the ray and the z-axis. By comparing the magnitudes of θ, φ, and γ, the coordinate axis corresponding to the smallest angle is selected as the principal direction of the ray.

[0124] The specific process of applying the continuous medium approximation to the principal direction of light is as follows:

[0125] Along the principal direction of the ray, select one point each on the outer boundary, the shock wave, and the inner boundary of the shock domain. Construct a polynomial equation about the coordinates of the principal direction of the ray based on the coordinates of these three points and their refractive index values. Taking the y-axis as the principal direction, the distribution of the refractive index values ​​of the three points along the principal direction of the ray is as follows: , , Based on this, the analytical expression for the continuous refractive index distribution in the principal direction of the shock wave domain can be derived. Substitute the three points into... The polynomial equation is obtained:

[0126] ;

[0127] Solve the polynomial using Cramer's rule:

[0128] ;

[0129] ;

[0130] ;

[0131] ;

[0132] in, is the determinant of the coefficient matrix of the polynomial equation; For the basic determinant Based on this, The corresponding coefficient column is replaced with the value column of refractive index n. n 1, n 2, n 3) The determinant of the new matrix formed after that; For the basic determinant Based on this, The corresponding coefficient column is replaced with the value column of refractive index n. n 1, n 2, n 3) The determinant of the new matrix formed after that; For the basic determinant Based on this, The corresponding coefficient column is replaced with the value column of refractive index n. n 1, n 2, n3) after, the determinant of the new matrix formed;

[0133] Finally, the polynomial coefficients are solved 、 、 :

[0134] ;

[0135] ;

[0136] ;

[0137] At this time, the flow field in the main direction of the light ray is a continuous medium, and the flow field in the non-main direction of the light ray is in a discrete state.

[0138] The specific process of implementing the structured grid discrete processing on the non-main direction of the light ray is as follows:

[0139] In a two-dimensional flow field, the non-main direction of the light ray is divided into equal-interval micro-intervals, and in a single micro-interval, the flow field is continuous in the main direction.

[0140] In a three-dimensional flow field, the non-main direction of the light ray is divided into micro-units of equal size, and in a micro-unit, the flow field is also continuous in the main direction of the light ray.

[0141] The application can greatly reduce the number of complex interpolation operations by adopting the mixed dimension numerical continuous and discrete strategy, making the main direction of the light ray continuous and the non-main direction of the light ray discrete, fundamentally reducing the calculation complexity, realizing the rapid modeling of the light field of the hypersonic flow field, and being especially suitable for dynamic process analysis.

[0142] The application can greatly reduce the demand for computer memory and computing resources on the premise of ensuring spatial resolution by focusing on the key area (shock wave domain) and adopting an efficient mixed dimension strategy, which makes it possible to realize the light field modeling of a high-resolution flow field of the order of ten million grids on a general computing platform, reduces the technical threshold and application cost, and has high engineering application value.

[0143] S3: Based on the analogy of Fermat's principle and the energy conservation condition of the Mopet principle, the transmission of the light ray in the refractive index field is analogized to the movement of the equivalent real particles in the refractive index potential field, so as to establish the dynamic relationship between the equivalent real particles and the refractive index potential field.

[0144] The specific process of establishing the dynamic relationship between the equivalent real particles and the refractive index potential field is as follows:

[0145] The variational form of Fermat's principle is:

[0146] ;

[0147] wherein, is the refractive index, is the light path;

[0148] The principle of Moppert under the condition of energy conservation is:

[0149] ;

[0150] wherein, is the optical momentum, is the generalized coordinate;

[0151] The Fermat principle and the principle of Moppert under the condition of energy conservation have the same form, by analogy between the two, the transmission of light in the refractive index field is analogous to the motion of the equivalent real particles in the refractive index potential field, at this time, the mass of the equivalent real particles is set to 1, then the dynamic relationship between the equivalent real particles and the refractive index potential field is:

[0152] ;

[0153] wherein, is the speed of the equivalent real particles.

[0154] The application is based on the Hamiltonian optics principle, the propagation behavior of light wave in the time-varying refractive index medium is equivalent to the dynamic process of particles with mass in the potential field, a rigorous dynamic theoretical basis is provided for the light field tracking, and the modeling accuracy of the shock light field is improved.

[0155] S4: based on the dynamic relationship between the equivalent real particles and the refractive index potential field, the straight line equation of the light and the spatial distribution analytical expression of the inner boundary of the shock domain are solved, the initial intersection point of the light and the shock domain is solved, the transmission path of the light in the shock is simulated, and the modeling of the shock light field in the hypersonic flow field is realized.

[0156] The straight line equation of the light and the spatial distribution analytical expression of the inner boundary of the shock domain are solved, the initial intersection point of the light equation and the shock domain is solved, and the refractive index at the initial intersection point is .

[0157] In the two-dimensional flow field, the straight line equation of the light can be converted into , and the spatial distribution analytical expression of the inner boundary of the shock domain can be converted into , and the initial intersection point of the straight line propagation of the light and the inner boundary of the shock domain can be obtained by combining the two:

[0158] .

[0159] The straight line equation of the light ray in the three-dimensional flow field can be converted to The analytical expression for the spatial distribution of the inner boundary of the shock domain can be converted to Both can be combined to obtain the initial intersection point of the light ray along the straight line and the inner boundary of the shock domain:

[0160] .

[0161] Taking the y-axis as the main direction of the light ray for exemplary illustration, Figure 6 shows the transmission path of the simulated light ray in the shock in a two-dimensional flow field, Figure 7 shows the transmission path of the simulated light ray in the shock in a three-dimensional flow field.

[0162] As shown in Figure 6 , in a two-dimensional flow field, based on the dynamic relationship between the equivalent real particles and the refractive index potential field, the velocity component expressions of the equivalent real particles in the main direction and the non-main direction of the light ray are obtained:

[0163] ;

[0164] ;

[0165] wherein, is the velocity component of the equivalent real particle in the non-main direction of the light ray, is the velocity component of the equivalent real particle in the main direction of the light ray, is the refractive index distribution in the main direction of the light ray in a two-dimensional flow field, is the angle between the light ray and the x-axis at the initial intersection point.

[0166] The velocity component expressions of the equivalent real particles in the main direction and the non-main direction of the light ray derive the position differential expression of the equivalent real particles:

[0167] .

[0168] Integrating both sides of the position differential expression of the equivalent real particles respectively obtains the position expression:

[0169] .

[0170] When the equivalent real particle passes through a small interval in the non-main direction of the light ray to the next small interval, the refractive index n in the small interval is updated, and the motion trajectory of the equivalent real particle in the small interval is calculated based on the updated refractive index n.

[0171] As shown in Figure 7As shown, in the three-dimensional flow field, the velocity component expressions of the equivalent real particles in the non-principal direction plane and the principal direction of the light are obtained based on the dynamic relationship between the equivalent real particles and the refractive index potential field:

[0172]

[0173]

[0174] wherein, is the velocity component of the equivalent real particle in the x-z plane, is the velocity component of the equivalent real particle in the principal direction of the light, is the refractive index distribution along the principal direction of the light in the three-dimensional flow field, is the angle between the light and the x-z plane at the initial intersection point;

[0175] The position differential expression of the equivalent real particle is derived based on the velocity component expressions of the equivalent real particle in the x-z plane and the y axis:

[0176]

[0177] wherein, is the projection of the equivalent real particle in the x-z plane;

[0178] The position expression is obtained by integrating the position differential expression of the equivalent real particle on both sides:

[0179]

[0180] When the equivalent real particle passes through a micro-unit in the non-principal direction of the light to reach the next micro-unit, the refractive index n in the micro-unit is updated, and the motion trajectory of the equivalent real particle in the micro-unit is calculated based on the updated refractive index n.

[0181] The ray initial intersection point solving and transmission path simulation operations are repeatedly performed to obtain the transmission paths of multiple rays in the shock wave domain, and the shock wave light field modeling in the hypersonic flow field is realized based on the transmission paths of the multiple rays.

[0182] The optical path length of a single ray The optical path length of a single ray is calculated by the following formula:

[0183]

[0184] wherein, s i represents the propagation distance of the light in each tracking step, n i represents the refractive index of each tracking step. ​​​​

[0185] Optical path difference The expression is:

[0186] ;

[0187] in, Indicates all The average value.

[0188] wavefront aberration It can be expressed by the following formula:

[0189] ;

[0190] in, k = 2π / λ , λ wavelength of light 。

[0191] wavefront aberration It is to calculate the actual pupil function. The basis:

[0192] ;

[0193] in, Represents the ideal pupil function. This represents the phase perturbation factor. Indicates wavefront aberration. It represents the imaginary unit.

[0194] Complex amplitude function It is the actual pupil function Fourier transform:

[0195] .

[0196] Point spread function It is a complex amplitude function The square of the modulus:

[0197] ;

[0198] in, yes .

[0199] Optical transfer function It is a point spread function Fourier transform:

[0200] .

[0201] Modulation transfer function It is the optical transfer function The model:

[0202] .

[0203] So far, the light field modeling of the shock wave in the hypersonic flow field has been completed, and the numerical model of the interaction between the shock wave and the light field is constructed. On this basis, the optical transmission characteristics of the shock wave region are also calculated.

[0204] The high efficiency of the modeling method of the application enables it to seamlessly connect the nanosecond time step data generated by the computational fluid dynamics simulation, thereby realizing effective simulation and analysis of the dynamic and high-frequency evolution process of the aerodynamic optical effects (such as image blur and jitter) of the hypersonic flow field, and providing a powerful numerical simulation tool for the performance evaluation of the optical guidance and observation system of the aircraft.

[0205] It should be understood that the steps can be reordered, added or deleted using the various forms of the flow shown above. For example, the steps described in the present disclosure can be executed in parallel, sequentially or in a different order, as long as the desired results of the technical solutions of the present disclosure can be achieved, which is not limited herein.

[0206] The above specific embodiments do not constitute a limitation on the protection scope of the present application. Those skilled in the art should understand that various modifications, combinations, sub-combinations and substitutions can be made according to design requirements and other factors. Any modifications, equivalent replacements and improvements made within the spirit and principles of the present application shall be included in the protection scope of the present application.

Claims

1. A method for modeling the shock wave optical field in a hypersonic flow field, characterized in that, Comprise the following steps: S1: using unsupervised learning algorithm to flow field data gradient clustering, determine the position of the shock wave, the inside and outside boundary of the shock wave field, and obtain the inside and outside boundary of the shock wave field, the spatial distribution of the shock wave is solved by numerical fitting method; S2: construct mixed dimension numerical continuous and discrete strategy, determine the main direction of the light based on the strategy, implement continuous medium approximation processing to the main direction of the light, implement structured grid discrete processing to the non-main direction of the light; S3: based on the analogy of Fermat principle and the analogy of the principle of momentum under the condition of energy conservation, the transmission of light in the refractive index field is analogous to the motion of equivalent real particles in the refractive index potential field, thereby establishing the dynamic relationship between the equivalent real particles and the refractive index potential field; In step S3, the specific process of establishing the dynamic relationship between the equivalent real particles and the refractive index potential field is: The variational form of Fermat principle is: ； wherein n is the refractive index, is the light path; The principle of momentum under the condition of energy conservation is: ； wherein is the optical momentum, is the generalized coordinate; By analogy of Fermat principle and the principle of momentum under the condition of energy conservation, the transmission of light in the refractive index field is analogous to the motion of equivalent real particles in the refractive index potential field, at this time the mass of the equivalent real particle is set to 1, then the dynamic relationship between the equivalent real particle and the refractive index potential field is established: ； wherein is the magnitude of the velocity of the equivalent real particle; S4: based on the dynamic relationship between the equivalent real particle and the refractive index potential field, the straight line equation of the light and the spatial distribution of the inside boundary of the shock wave field are solved, the initial intersection of the light and the shock wave field is solved, the transmission path of the light in the shock wave is simulated, and the modeling of the shock wave light field in the hypersonic flow field is realized.

2. The method of claim 1, wherein, In step S1, the specific process of using unsupervised learning algorithm to flow field data gradient clustering is: S11: abstract and simplify the flow field data into one-dimensional data set; wherein, the flow field data is composed of node coordinates and flow field attribute data, and the node position is constant; S12: randomly select two nodes in the flow field data as the initial centroid, calculate the density value of each node in the flow field data and the Euclidean distance between the two initial centroids, and classify the density value of each node according to the size of the Euclidean distance to the cluster corresponding to the nearest initial centroid; S13: recalculate the new centroid of each cluster to replace the old centroid, repeat steps S11 and S12 until the centroid and cluster assignment converge and no longer change, and divide the flow field into far field and near field; S14: take the demarcation of far field and near field as the outer boundary of shock wave field, perform downward translation operation on the outer boundary of shock wave field, when the refractive index gradient reaches the peak value, take the peak position as the position of shock wave; S15: select the node symmetrical to the position of shock wave in the outer boundary of shock wave field as the inner boundary of shock wave field.

3. The method of claim 1, wherein, In step S2, the determination method of the main direction of the light is: In two-dimensional flow field, compare the size of the angle between the light and x-axis, y-axis, select the coordinate axis with the smallest angle as the main direction of the light; In three-dimensional flow field, compare the size of the angle between the light and x-axis, y-axis, z-axis, select the coordinate axis with the smallest angle as the main direction of the light.

4. The method of claim 1, wherein, In step S2, the specific process of implementing continuous medium approximation processing to the main direction of the light is: A point is selected on the outer boundary of the shock domain, the shock, and the inner boundary along the main direction of the light ray, and a polynomial equation about the main direction coordinate is constructed according to the coordinates and the refractive index values of the three points. The polynomial equation is solved by the Cramer rule to obtain an analytical expression of the continuous distribution of the refractive index of the shock domain in the main direction of the light ray, and the continuous medium approximation processing of the main direction of the light ray is realized.

5. The method of claim 1, wherein, In step S2, the specific process of implementing the structured grid discrete processing on the non-main direction of the light ray is as follows: In a two-dimensional flow field, the non-main direction of the light ray is divided into equal-interval micro-intervals, and the flow field in a single micro-interval remains continuous in the main direction of the light ray; In a three-dimensional flow field, the non-main direction of the light ray is divided into micro-units of equal size, and the flow field in a single micro-unit remains continuous in the main direction of the light ray.

6. The method of claim 1, wherein, In step S4, let the refractive index of the initial intersection be The specific process of the transmission path of the simulation light ray in the shock domain is as follows: If the flow field is a two-dimensional flow field and the main direction of the light ray is the y-axis, the velocity component expressions of the equivalent physical particles in the main direction and the non-main direction of the light ray are obtained based on the dynamic relationship between the equivalent physical particles and the refractive index potential field: ； ； wherein is the velocity component of the equivalent real particle in the non-principal direction of the light ray, is the velocity component of the equivalent real particle in the principal direction of the light ray, is the refractive index distribution in the two-dimensional flow field along the principal direction of the light ray, is the angle of the light ray with the x-axis at the initial intersection point; The position differential expression of the equivalent physical particles is derived from the velocity component expressions of the equivalent physical particles in the main direction and the non-main direction of the light ray: ； The position expression of the equivalent physical particles is obtained by integrating the two sides of the position differential expression of the equivalent physical particles: ； When the equivalent physical particle passes through a micro-interval in the non-main direction of the light ray, the refractive index n in the micro-interval is updated, and the motion trajectory of the equivalent physical particle in the micro-interval is calculated based on the updated refractive index n; If the flow field is a three-dimensional flow field and the main direction of the light ray is the y-axis, the velocity component expressions of the equivalent physical particles in the main direction and the non-main direction of the light ray are obtained based on the dynamic relationship between the equivalent physical particles and the refractive index potential field: ； ； wherein is the velocity component of the equivalent real particle in the x-z plane, is the velocity component of the equivalent real particle in the main direction of the light ray, is the refractive index distribution along the main direction of the light ray in the three-dimensional flow field, is the angle of the light ray with the x-z plane at the initial intersection point; The position differential expression of the equivalent physical particles is derived from the velocity component expressions of the equivalent physical particles in the x-z plane and the y-axis: ； wherein is the projection of the equivalent real particle on the x-z plane; The position expression of the equivalent physical particles is obtained by integrating the two sides of the position differential expression of the equivalent physical particles: ； When the equivalent physical particle passes through a micro-unit in the non-main direction of the light ray, the refractive index n in the micro-unit is updated, and the motion trajectory of the equivalent physical particle in the micro-unit is calculated based on the updated refractive index n.

7. The method of claim 1, wherein, In step S4, the specific process of solving the initial intersection point of the light ray and the shock domain by combining the straight line equation of the light ray with the analytical expression of the spatial distribution of the inner boundary of the shock domain is as follows: In the two-dimensional flow field, the straight line equation of the light ray is converted to The spatial distribution analytical expression of the inner boundary of the shock wave domain is converted to The initial intersection point of the light ray along the straight line and the inner boundary of the shock wave domain is obtained by combining the two. ； In the three-dimensional flow field, the straight line equation of the light ray is converted to The spatial distribution analytical expression of the inner boundary of the shock wave domain is converted to Both together can get the initial intersection point of the light ray along the straight line and the inner boundary of the shock wave domain: 。 8. The method of claim 1, wherein, In step S4, the specific process of realizing the modeling of the shock light field in the hypersonic flow field is as follows: The optical path length of each light ray is calculated The calculation formula is: ； wherein s i is the propagation distance of the light ray in each tracing step, n i is the refractive index of each tracing step; calculating the optical path difference of each light ray based on the optical path length of each light ray , and the calculation formula is: ； wherein the average of all the average of all based on the optical path difference of each light ray , the wavefront aberration is calculated , the calculation formula is: ； wherein k=2π / λ , λ is the wavelength of the light. Wavefront aberration based , the actual pupil function is calculated , the formula is: ； wherein is an ideal pupil function, denotes a phase perturbation factor; Based on the actual pupil function , the complex amplitude function is calculated as ； Based on complex amplitude function , the point spread function is calculated , the calculation formula is: ； wherein is the complex conjugate of Point spread function based Optical transfer function The formula is: ； Based on the optical transfer function , the modulation transfer function , the formula is: 。 9. The method of claim 8, wherein, The transmission paths of multiple light rays in the shock are simulated, and the modeling of the shock light field in the hypersonic flow field is completed based on the transmission paths of the multiple light rays.

Images