Nothing Special   »   [go: up one dir, main page]

Next Article in Journal
Novel Corrugated Long Period Grating Surface Balloon-Shaped Heterocore-Structured Plastic Optical Fibre Sensor for Microalgal Bioethanol Production
Previous Article in Journal
End-Point Position Estimation of a Soft Continuum Manipulator Using Embedded Linear Magnetic Encoders
You seem to have javascript disabled. Please note that many of the page functionalities won't work as expected without javascript enabled.
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Quick Simulation Method for Aero-Optical Effects Based on a Density Proxy Model

1
School of Astronautics, Beihang University, Beijng 100191, China
2
Beijing Institute of Control and Electronic Technology, Beijing 100038, China
3
Qian Xuesen Laboratory of Space Technology, Beijing 100094, China
*
Author to whom correspondence should be addressed.
Sensors 2023, 23(3), 1646; https://doi.org/10.3390/s23031646
Submission received: 10 December 2022 / Revised: 21 January 2023 / Accepted: 30 January 2023 / Published: 2 February 2023
(This article belongs to the Section Optical Sensors)
Figure 1
<p>Schematic diagram of boundary layer density proxy model.</p> ">
Figure 2
<p>Schematic diagram of simulation threshold of density proxy model.</p> ">
Figure 3
<p>Schematic diagram of ellipsoidal vortex coordinate system.</p> ">
Figure 4
<p><span class="html-italic">SR</span> of density proxy model changing with control parameters.</p> ">
Figure 5
<p>The solution results of the hyper-parameters based on Bayesian optimization: (<b>a</b>) iterative scatter distribution of <span class="html-italic">N</span><sub>1</sub>, <span class="html-italic">N</span><sub>2</sub>; (<b>b</b>) iterative scatter distribution of <span class="html-italic">N</span><sub>3</sub>, <span class="html-italic">N</span><sub>4</sub>; (<b>c</b>) iteration trajectory of the trust coefficient <math display="inline"><semantics> <mrow> <msub> <mi>γ</mi> <mi>R</mi> </msub> </mrow> </semantics></math>; (<b>d</b>) the change curve of loss function.</p> ">
Figure 6
<p>(<b>a</b>) <math display="inline"><semantics> <mrow> <mi>S</mi> <msup> <mi>R</mi> <mi>s</mi> </msup> </mrow> </semantics></math> of the turbulence density proxy model; (<b>b</b>) <math display="inline"><semantics> <mrow> <mi>S</mi> <msup> <mi>R</mi> <mi>p</mi> </msup> </mrow> </semantics></math> of the distortion prediction model.</p> ">
Figure 7
<p>Optical characteristics of XY plane in density proxy model: (<b>a</b>) <math display="inline"><semantics> <mrow> <mi>t</mi> <mo>=</mo> <mn>0</mn> <mo> </mo> <mi>ms</mi> </mrow> </semantics></math>; (<b>b</b>) <math display="inline"><semantics> <mrow> <mi>t</mi> <mo>=</mo> <mn>5</mn> <mo> </mo> <mi>ms</mi> </mrow> </semantics></math>; (<b>c</b>) <math display="inline"><semantics> <mrow> <mi>t</mi> <mo>=</mo> <mn>10</mn> <mo> </mo> <mi>ms</mi> </mrow> </semantics></math>; (<b>d</b>) <math display="inline"><semantics> <mrow> <mi>t</mi> <mo>=</mo> <mn>15</mn> <mo> </mo> <mi>ms</mi> </mrow> </semantics></math>.</p> ">
Figure 8
<p>(<b>a</b>) OPD of the turbulence density proxy model; (<b>b</b>) OPD of the flow field by the CFD.</p> ">
Figure 9
<p>(<b>a</b>) OPD spectrum distribution of flow direction; (<b>b</b>) OPD spectrum distribution error of flow direction.</p> ">
Figure 10
<p>Monte Carlo verification results under different conditions: (<b>a</b>) distribution of sampling points at 10 km; (<b>b</b>) relative error of OPD at 10 km; (<b>c</b>) distribution of sampling points at 15 km; (<b>d</b>) relative error of OPD at 15 km.</p> ">
Figure 11
<p>Simulation of distorted star maps: (<b>a</b>) star A at 15 km; (<b>b</b>) star B at 15 km; (<b>c</b>) star A at 5 km; (<b>d</b>) star B at 5 km.</p> ">
Figure 12
<p>The simulation results of distorted star maps with the previous study [<a href="#B27-sensors-23-01646" class="html-bibr">27</a>].</p> ">
Review Reports Versions Notes

Abstract

:
Aero-optical effects caused by high-speed flow fields will interfere with the transmission of starlight, reduce the accuracy of optical sensors, and affect the application of celestial navigation on hypersonic vehicles. At present, the research of aero-optical effects relies heavily on the flow field simulation of computational fluid dynamics (CFD), which requires a great deal of computing resources and time, and cannot satisfy the demand of the rapid analysis of aero-optical effects in the engineering design stage. Therefore, a quick simulation method for aero-optical effects based on a density proxy model (DP-AOQS) is proposed in this paper. A proxy model of the turbulent density field is designed to replace the density field in the CFD simulation, and the proxy model is parametrically calibrated to simulate the optical characteristics of the turbulent boundary layer (TBL) in the external flow field of the optical window. The performance of DP-AOQS in the visible light band is verified from the perspectives of density field distribution, optical path difference (OPD), and fuzzy star map. The simulation results show that the method can quickly provide the distortion results of aero-optical effects in different flight conditions on the premise of ensuring the simulation accuracy. The research in this paper provides a new analytical method for the study of aero-optical effects.

1. Introduction

With the rapid development of the aerospace industry in various countries, hypersonic vehicles have become a major research focus of aircraft technology [1,2]. In order to improve the autonomy and accuracy of navigation and guidance, hypersonic vehicles are often equipped with optical imaging detection equipment to provide detection information for navigation and guidance. The application of star sensors in celestial navigation is a particularly extensive research direction at present [3,4]. However, when hypersonic vehicles fly at a high speed in the atmosphere, complex flow structures, such as turbulent boundary layers, shock waves, and expansion waves and mixed layers, may be formed outside the optical window, which makes the density field, temperature field, and velocity field of the air change dramatically [5,6]. When light passes through the complex flow field, the high-speed flow field will disturb the transmission of light and cause optical distortion; that is, aero-optical effects which can make the target image formed by imaging detection equipment appear subject to energy attenuation, offset, blur, jitter, and other phenomena [7,8,9,10]. Therefore, the research of aero-optical effects is of great significance for the application of optical sensors in hypersonic vehicles.
The law of light disturbance caused by high-speed flow field is important to lay a foundation for schemes designed to reduce the influence of aero-optical effects. At present, the main method of research is to obtain the external flow field data of hypersonic vehicles by using high-precision numerical calculation methods, such as computational fluid dynamics (CFD) and direct simulation Monte Carlo (DSMC) [11,12,13]. Then, the optical distortion and image quality evaluation of optical imaging equipment can be obtained through geometric optics, physical optics, statistical optics, and other optical transmission calculation methods [14,15,16]. The focus of those methods is mainly on how to obtain the flow field data, which requires a huge amount of computation and time. When analyzing the influence of flight conditions on the aero-optical effects, most of them adopt the method of repeatedly solving the flow field to calculate the aero-optical effects under the corresponding conditions [17,18,19]. Some scholars have studied the vortex structure of the boundary layer flow field by the CFD method, but the solution of complex differential equation consumes a lot of computational resources [20]. Although the CFD and DSMC methods have a high accuracy in calculating the density field, and the aero-optical effect results calculated from them have a high reliability, the calculation takes a long time and the utilization rate of the data is low. Due to the complexity of turbulence, it is necessary to recalculate the full flow field data when the flight conditions change slightly. Numerical simulation methods such as these, with a high precision and high computation, are used to calculate the density field to analyze the aero-optical effects, which cannot satisfy the requirement of rapidity when conducting aero-optical engineering research at the aircraft design stage. Although some scholars have improved various numerical methods based on hydrodynamics [21], in terms of the aero-optical effects, we mainly focus on the optical properties of the flow field. Therefore, it is necessary to find a quick simulation method for the aero-optical effects expected for the research method based on CFD.
As the calculation of a complex flow field is the main factor affecting computation in aero-optical simulations, a proxy model of the flow field can effectively improve the simulation efficiency of aero-optical effects. Many scholars have verified, through many experiments and numerical simulations, that the aero-optical effects have a corresponding relationship with the boundary layer thickness, Mach number, incoming flow density and other parameters of the optical window boundary layer, and obtained scaling law models such as the root mean square of the optical path difference and jitter angle [22,23,24,25,26]. These high-fidelity aero-optical distortion laws are the basis for matching the flow-field proxy model with the actual high-speed flow field.
The contribution of this paper is to design a boundary layer density proxy model to replace the density field in CFD numerical simulation and realize the rapid simulation of aero-optical effects, as shown in Figure 1. In this paper, the density proxy model is based on the density ellipsoidal vortex model in the aero-optical effects steady-state simulator proposed previously [27]. In contrast to the previous method, this paper considers the density distribution inside the spherical vortex and converts the large-scale vortex structure into the superposition of multiple ellipsoidal vortices. To conform to the real high-speed flow field environment, the motion model of the ellipsoidal vortex structure is added on the basis of the previous static structure to make the density proxy model closer to the real high-speed flow field. In order to ensure the accuracy of the model, polynomial and Bayesian optimization methods are used to calibrate the model. The density and scale control parameters are calibrated based on the boundary layer distortion prediction model. The disturbance parameters solved by the boundary layer linear stability theory are used as the motion parameters of the proxy model. The control parameters of the flow direction movement and dip distribution are calibrated based on the turbulent density fluctuation distribution law. Finally, a density proxy model reflecting the optical characteristics of high-speed flow field is obtained, and the fast simulation of aero-optical effects is completed through the calculation of light.
The remainder of this paper is organized as follows: in Section 2, the design of the density proxy model in DP-AOQS is presented. The calibration method of the control parameters in the density proxy model is described in Section 3. In Section 4, the aero-optical effects in the visible light band are verified from the perspectives of density field distribution, OPD, and fuzzy star map. Finally, the conclusions are drawn in Section 5.

2. Design of Density Proxy Model

In this section, the density proxy model of the boundary layer flow field is established. This paper mainly studies the boundary layer flow field outside the optical window, and the thickness of the TBL under hypersonic conditions is about 5–20 mm. In addition, considering the size of the optical window, in order to fully simulate the transmission disturbance of the incident starlight by the TBL flow field outside the optical window fully, a cuboid space above the optical window was selected as the simulation domain with the size of L x × L y × L z = 200   mm × 30   mm × 100   mm , as shown in Figure 2. OXYZ is defined as the window coordinate system. The coordinate origin O is located at the centroid of the upper surface of the optical window.
The density sphere vortex model was used to establish the turbulence density proxy model, in which the density sphere vortex model is used to simulate the large-scale structure in turbulence that is much larger than the wavelength of the incident light and has density spatial coherence.

2.1. Continuous Density Model of Ellipsoidal Vortex

The density ρ V ( r , t ) of the large-scale vortex structure is composed of multiple ellipsoidal vortex densities, which can be expressed as follows:
ρ V ( r , t ) = i ρ i ( r , t )
where ρ i ( r , t ) is the density distribution of the i -th ellipsoidal vortex, which is different from the uniform ellipsoidal vortex in the literature [27]. In order to ensure the continuity of the density field, the internal density ρ i ( r i , ψ i , h i ) of the ellipsoidal vortex model is designed as a layered distribution for time t :
ρ i ( r i , ψ i , h i ) = ρ i n × r i Λ h i / Λ h i 2 h i 2
where r i , ψ i and h i are the cylindrical coordinates under the ellipsoidal vortex coordinate system o i x i , V y i , V z i , V , as shown in Figure 3.
o i is the centroid of the i -th ellipsoidal vortex; axis x i , V points to the major axis of the ellipsoid and is located in plane OXY; axis o i y i , V points to the minor axis; and Λ h i is the major axis length of the ellipsoidal surface at a distance h i from the plane o i x i , V y i , V . ρ in is the maximum spatial fluctuation in the density inside the ellipsoidal vortex, which can be expressed as follows:
ρ in = 1 N k ρ F M a , y Λ ρ + ρ
where ρ is the incoming flow density under the current conditions. 1 N represents that the density change inside the ellipsoidal vortex may be higher or lower than the incoming flow density. Λ is the major axis length of the ellipsoidal vortex and k ρ is the control parameter of the ellipsoidal vortex density. F M a , y is the Mach number influence factor, which is given by:
F M a , y = f y / δ M a 2 1 + γ 1 2 M a 2 1 / 2 1 + γ 1 2 1 f 2 y / δ M a 2 3 / 2
where f y / δ = U y / δ / U is the average flow velocity distribution; U y / δ is the flow velocity at the height y / δ , which can be approximated by the logarithmic law c 1 ln y / δ + c 2 in a completely turbulent state; and U is the incoming flow velocity.

2.2. Location and Scale Model of Ellipsoidal Vortex

The position and scale distribution of the ellipsoidal vortex structure in the density proxy model are modeled by using the increasing function relationship. The position distribution of the ellipsoidal vortex is established as a Gaussian distribution model about the wall height:
P y _ = K exp y _ y _ 0 2 2 σ 2
where y _ = y / δ and P y _ are the probability of arranging an ellipsoidal vortex at the wall height y . y _ 0 is the wall height with the maximum probability of placing the ellipsoidal vortex. Results from different studies have shown that structures at a height of approximately 0 . 8 δ contribute most to optical distortions [14,28], so we took y _ 0 = 0.8 δ . σ is the variance of Gaussian distribution, and we took it as a fixed value of 0 . 5 δ and used the size of the gas ellipsoid to adjust the shape of the distorted wavefront. The normalization coefficient K satisfies:
K i = 1 k n u m exp y _ y _ 0 2 2 σ 2 = 1
where k n u m is the total number of ellipsoidal vortices in the simulation domain, and the typical value is between hundreds and thousands [29], so k n u m = 5000 was taken in this paper. The size of the large-scale structure in the boundary layer increases at a certain height. The major axis length of the ellipsoidal vortex is described by a quadratic function as follows:
Λ y _ = k l e n y _ y _ min 2 + Λ min
where y _ min is the wall height of the smallest ellipsoidal vortex; namely, the height for the vertex of the quadratic function, taken as 0.05 δ according to the literature [14]. Λ min is the minimum value of Λ y _ . The range of the minimum density correlation length in TBL is approximately 0.05 δ ~ 0.1 δ [14,28], so we took Λ min = 0.05 δ . k l e n is taken as the control parameter of the gas-ellipsoidal scale. The density correlation length range of the turbulent boundary layer is about 0.15 δ Λ ( y _ ) max 0.4 δ [17,27]. The constraint range of k l e n can be obtained according to Equation (7).

2.3. Motion Model of Ellipsoidal Vortex

There are disturbance waves with different wave numbers and phase velocities in the supersonic boundary layer [30]. The motion of the ellipsoidal vortex model is modeled as disturbance waves, in which the centroid ( x i , y i , z i ) of the ellipsoidal vortex moves in the form of the superposition of the flow direction traveling waves:
x ˙ i = k u ^ k exp j ( α k x i + β k z i ω k t + φ k ) + c . c
where u ^ k is the traveling wave disturbance characteristic function; α k is the k -th order disturbance flow beam; β k is the k -th order disturbance spread beam; ω k is the k -th order disturbance frequency; φ k is the k -th order mode phase; and j is an imaginary unit. For the time mode, α k is a real number and ω k is a complex number; that is, ω k = ω k , r + ω k , j j . c . c is a complex conjugate.
In the hypersonic boundary layer, the first mode has an important influence on the transition, and the second mode disturbance becomes the most unstable disturbance wave after the incoming Mach number exceeds four. Therefore, the ellipsoidal vortex center motion is set as the superposition of the first mode and the second mode; that is, k = 2 . The first mode flow direction beam α 1 is 0.10892184, the spread direction beam β 1 is ± 0.313 , the second mode flow direction beam α 2 is 1.44368594, and β 2 is ± 1.5 . According to the linear stability theory, the disturbance characteristic functions u ^ k and ω k in the control parameters of the ellipsoidal vortex motion can be obtained by solving the Orr–Sommerfeld small-disturbance equation. Considering the directional dependence of the large-scale turbulent structure, the expected inclination of the ellipsoidal vortex was set to 24.17°, and the oblateness was set to 0.447.

3. Calibration of Control Parameters in Density Proxy Model

In this section, the control parameters of the density proxy model are calibrated. k ρ , k l e n was used to control the position distribution of the ellipsoidal vortex and the spatial fluctuation of the density in the density proxy model, which can be calibrated by the methods of polynomial and Bayesian optimization.

3.1. Density and Scale Control Parameter Constraints

The density control parameter k ρ affects the density spatial fluctuation of the ellipsoidal vortex structure, and the scale control parameter k l e n affects the size distribution of the ellipsoidal vortex, both of which are closely related to the incoming flow state. The density correlation length range of the turbulent boundary layer is about 0.15 δ ~ 0.4 δ [17,27]. The constraint range of k l e n can be obtained according to Equation (7). To determine the constraint range of k ρ , the Strehl Ratio (SR) of the density proxy model was calculated in a steady state for different control parameters under the condition of a flight altitude of 10 km and Mach number of 5 Ma, as shown in Figure 4.
It can be seen from Figure 4 that S R s gradually decreases from 1 to 0 with the increase in the control parameters k ρ and k l e n . When k ρ 0.001 , the change in S R s with the increase in k l e n can be ignored, and the turbulence density proxy model cannot reflect the spatial fluctuation of turbulence. When k ρ 100 , S R s decreases rapidly with the increase in k l e n and the decreasing speed continues to increase with the increase in k ρ , resulting in the excessive spatial pulsation of the turbulence proxy model. Therefore, the upper limit k ρ b and lower limit k ρ a are used for the density control parameter k ρ , meeting the following conditions:
k ρ a = arg max k ρ   S R s ( k ρ , k l e n b ) S R a k ρ b = arg min k ρ   S R s ( k ρ , k l e n a ) S R b , k ρ k ρ a
where S R a and S R b are the minimum and maximum values of SR under the working conditions to be calibrated, respectively.

3.2. Calibration of Density and Scale Control Parameters

Under the given inflow conditions, the mean square value of OPD of the TBL can be described by the boundary layer distortion prediction model [25]:
OPD r m s = B K G D ρ M a 2 δ C f F M a
where B is a constant, taken as 0.2, F M a is an empirical function of the incoming Mach number, and C f is the wall friction coefficient of the TBL of the flat plate, which is related to the Reynolds number R e x . The density fluctuation and scale structure in turbulence are affected by the Mach number and Reynolds number of the incoming flow. The density control parameters affect the density fluctuation of the ellipsoidal vortex structure, and the scale control parameters affect the size distribution of the ellipsoidal vortex, both of which determine the optical characteristics of the density proxy model. Therefore, the control parameters k ρ and k l e n can be calibrated as a function of M a , R e x , and flight altitude H , which can be expressed as follows:
k l e n = g l e n ( M a , R e x , H ) k ρ = g ρ ( M a , R e x , H )
The working condition is expressed as y = M a R e x H T . g l e n ( M a , R e x , H ) and g ρ ( M a , R e x , H ) are expressed as g ( y ) . With the boundary layer distortion prediction model as the constraint condition, the optimal function mapping g is obtained; that is, a variational problem is solved:
g = arg min g O P D r m s p ( y ) O P D r m s s ( g ( y ) ) 2 2
where O P D r m s p ( y ) is the root mean square of the OPD of the boundary layer distortion prediction model under the given working condition y , and O P D r m s s ( g ( y ) ) is the root mean square of the OPD of the density proxy model. Because the gradient information of the control parameter output by the density proxy model cannot be obtained, the traditional variational method cannot be used to solve the variational problem in Equation (12). We transformed it into a parameter optimization problem, which is solved by a gradient-free parameter optimization algorithm.
SR in a density proxy model and distortion prediction model can be described as the form of multiplication or the addition of two univariate continuous functions with variable separation. Considering the existence of the multiplication form of the Mach number and Reynolds number in the distortion prediction model, we conducted research in the form of multiplication:
S R = h 1 ( x 1 ) h 2 ( x 2 )
where x 1 and x 2 in the density proxy model are k l e n and k ρ , respectively, and h 1 ( x 1 ) and h 2 ( x 2 ) are the monotonic increasing functions. x 1 and x 2 in the distortion prediction model are M a and R e x , respectively, h 2 ( x 2 ) is the monotonic decreasing function, and the monotonicity of h 1 ( x 1 ) is uncertain. Therefore, the control parameters k l e n and k ρ can be calibrated as the multiplication of two continuous functions with one variable separated:
k l e n = g l e n , 1 ( M a ) g l e n , 2 ( R e x ) k ρ = g ρ , 1 ( M a ) g ρ , 2 ( R e x )
The univariate continuous functions g l e n , i and g ρ , i can be fitted by the power function polynomial. The polynomial fitting function is set as g ( x ; a , N ) , where N is the highest order term of the power function and a = a 0 a 1 a N T is the coefficient vector of the power function.
g ( x ; π , N ) = k = 0 N a k x k
The coefficient a 1 , a 2 , a 3 , a 4 of the power function polynomial and the highest degree term N 1 , N 2 , N 3 , N 4 of the power function polynomial are parameters to be calibrated. For the convenience of writing, A and N are represented as A = a 1 a 2 a 3 a 4 and N = N 1 N 2 N 3 N 4 T , respectively. Then, Equation (11) can be written as:
k l e n = g ( M a ; a 1 , N 1 ) g ( l o g ( R e x ) ; a 2 , N 2 ) k ρ = g ( M a ; a 3 , N 3 ) g ( l o g ( R e x ) ; a 4 , N 4 )
The surface shape of S R p distribution with M a and R e x in the distortion prediction model changes obviously with the height. When fitting the functional relationship between k ρ , k l e n and M a , R e x by polynomial, it is necessary to calibrate the power function coefficient matrix A and the highest-order term N by height H . In order to determine the range of coefficient A to be optimized, the input variable M a , R e x , and the output variable k ρ , k l e n are normalized, and the value range is limited to the value interval [ 0 , 1 ] . M a and l o g ( R e x ) are calibrated according to the range of working conditions. Considering the scope of application of the distortion prediction model, the value range of M a is limited to 1.5 , 5 ; R e x is limited to 10 6 , 10 8 ; k l e n is limited to k l e n a / γ R , k l e n b / γ R ; and k ρ is limited to k ρ a / γ R , k ρ b / γ R . γ R is a trust coefficient which is used to change the range of the control parameters to be calibrated. Additionally, γ R is considered as a hyper-parameter, and the value range is set to 0.5 , 2 . When γ R < 1 , the calibration range of the control parameters is reduced, and when γ R > 1 , the calibration range of the control parameters becomes larger.
After the normalization of the input variables M a , l o g ( R e x ) and output variables k l e n , k ρ , the value of the higher-order term a i j of the power function polynomial coefficient is naturally limited to 1 , 1 , and the range of the zero-order power function coefficient A is 0 , 1 . Thus, coefficient matrix A satisfies the following equation:
A a A A b A a = 0 0 0 0 0 1 0 1 1 1 1 1 1 1 1 1 , A b = 1 1 1 1 1 0 1 0 1 1 1 1 1 1 1 1
The minimum value of the highest power term of the constrained power function is 3 and the maximum value is 20. The parameters A and N are calibrated by iterative optimization. For the fixed hyper-parameters N and γ R at the given flight altitude H 0 , the parameter A is optimized and solved. With the boundary layer distortion prediction model as the constraint, the optimal solution A of the parameter A to be calibrated is an optimization problem:
A = arg min A M , R e x O P D r m s p ( M a , R e x ) H 0 O P D r m s s ( M a , R e x ; A ) H 0 2 2 s u b j e c t t o A a A A b
where O P D r m s s ( M a , R e x ; A ) is the root mean square of the simulated OPD of the density proxy model at flight altitude H 0 . For a given flight altitude H 0 , the parameter N and trust coefficient γ R can be regarded as the hyper-parameters for optimization, and the loss function is defined as:
L o s s = d e f min A M , R e x O P D r m s p ( M a , R e x ) H 0 O P D r m s s ( M a , R e x ; A ) H 0 2 2 + a P log ( n i t e r )
where n i t e r is the total number of iterations in the optimization, a P log ( n i t e r ) is the penalty function of the total number of iterations, and a P is the penalty coefficient, taken as 0.001. When the number of iterations is too large, it means that the iterative optimization effect of the polynomial parameters is poor, and the polynomial form is not suitable as a fitting function. For the above hyper-parameter optimization problem, the Bayesian optimization algorithm is used to solve it:
x * = arg max x Χ f ( x )
where x is the parameter to be optimized, and its value range is Χ . f ( x ) is the objective function, and its gradient information is unknown. The basic idea of Bayesian optimization is to estimate the posterior distribution of the objective function based on the Bayesian principle. It is usually necessary to make assumptions about the model of the posterior distribution, estimate the posterior distribution on the basis of the observed data, and select new sampling points by maximizing the effect function to achieve the posterior distribution approaching the true distribution [31].
For H 0 = 15   km , the solution results of the hyper-parameters based on Bayesian optimization are shown in Figure 5, where (a) and (b) are the iterative scatter distribution of the hyper-parameter N . The red point is the optimal value, and the blue point is the initial iteration point. The size of the point increases with the number of iterations. Figure 5c shows the iteration trajectory of the trust coefficient γ R , which reaches the optimal value after 30 iterations. The change curve of the loss function L o s s is shown in Figure 5d. The prismatic sign is the loss function observed at the sampling point in the current iteration. The circle sign is the minimum value of the loss function in the current iteration. When the number of iterations is 12, the minimum value of the loss function is close to convergence. When the minimum value of the loss function reaches the tolerance range in 30 iterations, the Bayesian optimization solution is completed.
S R s of the turbulence density proxy model and S R p of the distortion prediction model after calibration are shown in Figure 6a,b. After polynomial calibration, the fitting error S R s S R p of SR is within 2 × 10 3 .

4. Simulation and Analysis

In this section, with a star sensor as the object, the simulation performance of DP-AOQS is verified from the characteristics of the density field distribution, OPD, and fuzzy star map of the density proxy model.

4.1. Verification of Optical Characteristics of Density Proxy Model

For the simulation conditions of H = 15   km , M a = 3 , R e x = 10 6 , the density fluctuation, and OPD of the density proxy model in the XY plane after being calibrated in Section 3 are shown in Figure 7. It can be seen from Figure 7 that the OPD fluctuated greatly with the flow direction distance because there were ellipsoidal vortices with uneven density fluctuations in the flow direction. For example, at t = 5   ms , there were wave peaks in the OPD near x = 0.08   m , and there were many ellipsoidal vortices (near y = 0.02   m 0.03   m ) with positive density fluctuations near x = 0.08   m . With the passage of time, the movement of the vortex structure in the direction of flow made the distribution of OPD fluctuate in both the time and space dimensions, which is consistent with the pulsation characteristics of the turbulent boundary layer [32]. It only took 200 s for DP-AOQS to obtain a set of results (simulation environment: Matlab R2021b, Intel i9-9900X CPU, 128 G RAM). A high-precision flow field structure calculated by CFD required 10–20 h (simulation environment: a high-performance server with AMD EPYC ROME 7H12, 128-cores CPU, 256 G RAM, the number of grids was more than 6 million, and the generation of large eddy simulation (LES) was more than 20,000).
In order to compare with the OPD calculated by the CFD, the OPD and projection surface of the density proxy model were obtained under the condition H = 20   km , M a = 3.8 , as shown in Figure 8a. The OPD of the flow field near the concave window under the same simulation conditions by the CFD is shown in Figure 8b [19]. It can be seen that there were also many peaks and troughs in the optical path difference surface. Figure 8 demonstrates that DP-AOQS can reflect the variation trend of the optical properties of the actual turbulent boundary layer flow field.
The flow direction OPD caused by the density fluctuation of the turbulent structure has certain amplitude frequency characteristics. Through experiments and similarity analysis, the amplitude spectrum model of the flow direction OPD was determined as follows [25]:
OPD x S t δ x = θ ^ p e a k 1 1 + S t δ x / 0.75 5 / 3 M a 2 f s ρ ρ S L δ U 2
where S t δ x is the normalized frequency of flow direction, satisfying S t δ x = δ f s / U . f s is the sampling frequency, ρ S L is the reference density of the sea level, ρ and U are the incoming flow density and velocity, respectively, and θ ^ p e a k is the peak value of the amplitude. To eliminate the difference caused by working conditions, we normalized Equation (21) to obtain the flow direction normalized spectrum OPD ^ x S t δ x , as shown in the following equation:
OPD ^ x S t δ x = OPD x S t δ x max ( OPD x S t δ x ) = 1 1 + S t δ x / 0.75 5 / 3
The amplitude spectrum model, density proxy model, and CFD flow field in the literature [19] were used to calculate the OPD spectrum distribution, as shown in Figure 9. It can be seen from Figure 9 that the density proxy model in this paper had a downward trend as the flow direction frequency increased, which is consistent with the trend of the amplitude spectrum model and the OPD spectrum obtained from the CFD calculation. However, the amplitude of the OPD spectrum in this paper decreased rapidly in the interval S t δ x = 10 1 10 0 because the flow direction motion model only used the first mode and the second mode for superposition, while there are many higher-order mode disturbance waves in the boundary layer. It is normal for the OPD spectrum to show buffeting in the interval S t δ x = 10 0 10 1 , which is caused by spatial discrete sampling, and amplitude mutation will occur at uncovered frequency points.

4.2. Monte Carlo Simulation under Different Working Conditions

The OPD of the density proxy model at 10 km and 15 km was verified by the Monte Carlo method, where the condition was M a 2 , 5 , R e x 1 × 10 6 , 1 × 10 8 , and the flight altitude was H H 0 1 , H 0 + 1 km , where H 0 is 10 km or 15 km. The results of the Monte Carlo verification are shown in Figure 10, where there were 100 sampling points in total.
From the relative error distribution of OPD in Figure 10, the fitting error of the density proxy model was within 10% at 15 km, while the flight altitude was 10 km, and the fitting error was within 12%, satisfying the accuracy requirements of the simulation. When the altitude decreases, the fitting error increases. This is because when the height decreases, the surface of SR will bend more violently with { M a , Re x } , and the error of the control parameters at the uncalibrated height increases. In general, the calibration of the control parameters of the density proxy model at a given altitude within 1 km could satisfy the accuracy requirements of OPD, and the calibration accuracy was higher than the steady-state simulation method proposed in the literature [27].

4.3. Simulation of Distorted Star Maps

After the starlight passes through the high-speed flow field, the star map on the image plane of the star sensor will degenerate. Because the distance between the star and the aircraft observed by the star sensor is far greater than the distance between the Earth and the sun, the observed star is regarded as a point light source at infinity. The star point signal on the CCD image plane is concentrated in a small circular area, where the energy distribution can be assumed to be the Gaussian distribution. The light field distribution on the image plane is the gray distribution of the diffraction image of a single star, and the display gray level of the corresponding pixel is obtained after gray level aggregation.
Under the condition of M a = 3 , R e x = 1 × 10 6 , the simulation results of the two distorted stars at 15 km and 5 km were as shown in Figure 11. In the simulation, the focal length of the star sensor was f = 61.88   mm , the pupil size was 50 mm, the pixel size was 15μm, the field angle was F O V = 20 , and the number of pixels on the image plane was 1456 × 1456 .
The gray value at the center of the distorted star map is listed in Table 1. When the flight altitude was 15 km, the S R p of the distortion prediction model was 0.8613, and the S R s of the turbulence density proxy model was 0.8650. At this point, compared with the standard star map, the gray distribution of star A and star B in the distorted star map was still the center point, and the energy was mainly concentrated in the range of one pixel in each direction from the center point. The gray value of the center point in the distorted star map decreased by 5–9 gray values.
When the flight altitude was 5 km, the S R p of the distortion prediction model was 0.0621, and the S R s of the turbulence density proxy model was 0.0645. Compared with the standard star map, the brightest point is still the center point of stars A and B in the distorted star map, but the energy distribution has spread in all directions, spreading 1–2 pixels. The gray value of the center point in the distorted star maps is greatly reduced, and the gray value of star A is the most obvious, falling 79 gray values. The distorted star maps obtained by DP-AOQS can be used for a reference in the design stage of celestial navigation in hypersonic vehicles. In addition, in order to highlight the advantages of this study, the distorted star map results of this paper are compared with the previous research methods. The star map obtained from the previous research [27] is shown in Figure 12.
By comparing the star map results in Figure 12 with those in Figure 11c,d, it can be found that the star map results are in line with objective laws, and the diffusion effect of the star map is more obvious because of the density change in the interior of the ellipsoid vortex and the motion model of the ellipsoid vortex considered in this paper.

5. Conclusions

In this paper, a quick simulation method of the aero-optical effects based on the density proxy model (DP-AOQS) is proposed. The density proxy model, considering the optical characteristics of the turbulent steady state and motion state, is designed to replace the density field in CFD numerical simulation. Based on the aero-optical laws obtained from the experiments and simulation analysis, the density proxy model is calibrated by combining polynomial fitting and Bayesian optimization so that the proxy model can reflect the optical characteristics in the actual high-speed flow field. The performance of the aero-optical effects in the visible light band is verified from the perspectives of the density field distribution, OPD, and fuzzy star map. Compared with the traditional CFD simulation method, DP-AOQS can quickly provide the distortion results of aero-optical effects under different flight conditions. It provides a powerful tool for aero-optical effects analysis in the design stage of celestial navigation in hypersonic vehicles.

Author Contributions

Conceptualization, B.Y.; methodology, H.Y. and Z.F.; project administration, H.Y.; validation, H.Y.; visualization, C.L.; writing—original draft, H.Y.; writing—review and editing, X.W. and J.M. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Science and Technology on Space Intelligent Control Laboratory of China (No. ZDSYS-2018-03), the National Natural Science Foundation of China (No. 61973018), and the Civil Aerospace Technology Pre-Research Project of China (No. D040301).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data sharing not applicable.

Acknowledgments

We thank the anonymous reviewers for their valuable comments which significantly improved the paper.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

ρ V ( r , t ) Density of the large-scale vortex structure.
ρ i ( r i , ψ i , h i ) Internal density of the ellipsoidal vortex model.
r i , ψ i , h i Cylindrical coordinates under the vortex coordinate system.
ρ Incoming flow density under the current conditions.
ρ in Maximum spatial fluctuation of the density inside the ellipsoidal vortex.
Λ Major axis length of the ellipsoidal vortex.
k ρ Control parameter of the ellipsoidal vortex density.
k ρ a Lower limit of k ρ .
k ρ b Upper limit of k ρ .
F M a , y Mach number influence factor.
f y / δ Average flow velocity distribution.
U y / δ Flow velocity at the height y / δ .
U Incoming flow velocity.
P y _ Probability of arranging an ellipsoidal vortex at wall height y .
y _ 0 Wall height with the maximum probability of placing ellipsoidal vortex.
k n u m Total number of ellipsoidal vortices in the simulation domain.
y _ min Wall height of the smallest ellipsoidal vortex.
Λ y _ Major axis length of the ellipsoidal vortex.
Λ min Minimum value of Λ y _ .
k l e n Control parameter of the gas-ellipsoidal scale.
k l e n a Lower limit of k l e n .
k l e n b Upper limit of k l e n .
u ^ k Traveling wave disturbance characteristic function.
α k The k -th order disturbance flow beam.
β k The k -th order disturbance spread beam.
ω k The k -th order disturbance frequency.
φ k The k -th order mode phase.
j Imaginary unit.
c . c Complex conjugate.
S R s Strehl ratio (SR) of the density proxy model.
OPD r m s Mean square value of OPD.
C f Wall friction coefficient of the TBL of the flat plate.
R e x Reynolds number.
H Flight altitude.
A Power function coefficient matrix.
N The highest order term of the power function.
a Coefficient vector of the power function.
γ R Trust coefficient which is used to change the range of control parameters.
L o s s Loss function.
n i t e r Total number of iterations in the optimization.
a P log ( n i t e r ) Penalty function of the total number of iterations.
a P Penalty coefficient.

References

  1. Quagliarella, D.; Pezzella, G.; Pirozzoli, S. An aerothermodynamic design optimization framework for hypersonic vehicles. Aerosp. Sci. Technol. 2019, 84, 339–347. [Google Scholar]
  2. Gruhn, P.; Gülhan, A. Aerodynamic measurements of an air-breathing hypersonic vehicle at Mach 3.5 to 8. AIAA J. 2018, 56, 4282–4296. [Google Scholar] [CrossRef]
  3. Yu, Y.J.; Xu, J.F.; Xiong, Z. SINS/CNS nonlinear integrated navigation algorithm for hypersonic vehicle. Math. Probl. Eng. 2015, 2015, 903054. [Google Scholar] [CrossRef]
  4. Bingbing, G.; Wenmin, L.; Gaoge, H.; Zhong, Y.; Xinhe, Z. Mahalanobis distance-based fading cubature Kalman filter with augmented mechanism for hypersonic vehicle INS/CNS autonomous integration. Chin. J. Aeronaut. 2022, 35, 114–128. [Google Scholar]
  5. Sun, X.; Yang, X.; Liu, W. Validation method of aero-optical effect simulation for supersonic turbulent boundary layer. AIAA J. 2021, 59, 410–416. [Google Scholar] [CrossRef]
  6. Mathews, E.; Wang, K.; Wang, M.; Jumper, E.J. Turbulence scale effects and resolution requirements in aero-optics. Appl. Opt. 2021, 60, 4426–4433. [Google Scholar] [CrossRef]
  7. Jumper, E.; Gordeyev, S. Physics and measurement of aero-optical effects: Past and present. Annu. Rev. Fluid Mech. 2017, 49, 419–441. [Google Scholar] [CrossRef]
  8. Wang, K.; Wang, M. Computational analysis of aero-optical distortions by flow over a cylindrical turret. AIAA J. 2016, 54, 1461–1471. [Google Scholar] [CrossRef]
  9. Emelyanov, V.; Teterina, I.; Volkov, K.; Yakovchuk, M. Aero-optical effects in free and wall-bounded turbulent compressible flows. Acta Astronaut. 2018, 150, 144–152. [Google Scholar] [CrossRef]
  10. Ding, H.; Yi, S.; Ouyang, T.; Shi, Y.; He, L. Influence of turbulence structure with different scale on aero-optics induced by supersonic turbulent boundary layer. Optik 2020, 202, 163565. [Google Scholar] [CrossRef]
  11. Zhang, D.; Tan, J.; Yao, X. Direct numerical simulation of spatially developing highly compressible mixing layer: Structural evolution and turbulent statistics. Phys. Fluids 2019, 31, 036102. [Google Scholar] [CrossRef]
  12. Zhang, B.; Liu, H.; Jin, S. An asymptotic preserving Monte Carlo method for the multispecies Boltzmann equation. J. Comput. Phys. 2016, 305, 575–588. [Google Scholar] [CrossRef]
  13. Wu, M.; Martin, M. Direct numerical simulation of supersonic turbulent boundary layer over a compression ramp. AIAA J. 2007, 45, 879–889. [Google Scholar] [CrossRef]
  14. Wang, K.; Wang, M. Aero-optics of subsonic turbulent boundary layers. J. Fluid Mech. 2012, 696, 122–151. [Google Scholar] [CrossRef]
  15. Gordeyev, S.; Jumper, E. Fluid dynamics and aero-optics of turrets. Prog. Aerosp. Sci. 2010, 46, 388–400. [Google Scholar] [CrossRef]
  16. Rennie, R.; Nguyen, M.; Gordeyev, S.; Jumper, E.; Cain, A.; Hayden, T. Wave-front measurements of a supersonic boundary layer using laser induced breakdown. AIAA J. 2017, 55, 2349–2357. [Google Scholar] [CrossRef]
  17. Tromeur, E.; Garnier, E.; Sagaut, P. Large-eddy simulation of aero-optical effects in a spatially developing turbulent boundary layer. J. Turbul. 2006, 7, N1. [Google Scholar] [CrossRef]
  18. Saxton-Fox, T.; McKeon, B.J.; Gordeyev, S. Effect of coherent structures on aero-optic distortion in a turbulent boundary layer. AIAA J. 2019, 57, 2828–2839. [Google Scholar] [CrossRef]
  19. Yang, B.; Yu, H.; Liu, C.; Wei, X.; Fan, Z.; Miao, J. A New Method for Analyzing the Aero-Optical Effects of Hypersonic Vehicles Based on a Microscopic Mechanism. Aerospace 2022, 9, 618. [Google Scholar] [CrossRef]
  20. Pirozzoli, S.; Bernardini, M.; Grasso, F. Characterization of coherent vortical structures in a supersonic turbulent boundary layer. J. Fluid Mech. 2008, 613, 205–231. [Google Scholar] [CrossRef]
  21. De Vanna, F.; Baldan, G.; Picano, F.; Benini, E. Effect of convective schemes in wall-resolved and wall-modeled LES of compressible wall turbulence. Comput. Fluids 2023, 250, 105710. [Google Scholar] [CrossRef]
  22. Sutton, G.W. Aero-optical foundations and applications. AIAA J. 1985, 23, 1525–1537. [Google Scholar] [CrossRef]
  23. Sutton, G.W. Effect of turbulent fluctuations in an optically active fluid medium. AIAA J. 1969, 7, 1737–1743. [Google Scholar] [CrossRef]
  24. Spina, E.; Smits, A.; Robinson, S. The Physics of supersonic turbulent boundary layers. Annu. Rev. Fluid Mech. 1994, 26, 287–319. [Google Scholar] [CrossRef]
  25. Gordeyev, S.; Smith, A.; Cress, J.; Jumper, E. Experimental studies of aero-optical properties of subsonic turbulent boundary layers. J. Fluid Mech. 2014, 740, 214–253. [Google Scholar] [CrossRef]
  26. Wyckham, C.; Smits, A. Aero-optic distortion in transonic and hypersonic turbulent boundary layers. AIAA J. 2009, 47, 2158–2168. [Google Scholar] [CrossRef]
  27. Bo, Y.; Zichen, F.; He, Y. Aero-optical effects simulation technique for starlight transmission in boundary layer under high-speed conditions. Chin. J. Aeronaut. 2020, 33, 1929–1941. [Google Scholar]
  28. Gao, Q.; Yi, S.; Jiang, Z.; He, L.; Wang, X. Structure of the refractive index distribution of the supersonic turbulent boundary layer. Opt. Lasers Eng. 2013, 51, 1113–1119. [Google Scholar] [CrossRef]
  29. Trolinger, J.; Rose, W. Technique for simulating and evaluating aero-optical effects in optical systems. In Proceedings of the 42nd AIAA Aerospace Sciences Meeting and Exhibit, Reno, NV, USA, 5–8 January 2004; pp. 471–481. [Google Scholar]
  30. Mack, L. Boundary-Layer Linear Stability Theory; Jet Propulsion Lab: Pasadena, CA, USA, 1969. [Google Scholar]
  31. Hoffman, M.; Shahriari, B. Modular Mechanisms for Bayesian Optimization. In NIPS Workshop Bayesian Optim. 2014, pp. 1–5. Available online: https://bayesopt.github.io/papers/2014/paper8.pdf (accessed on 9 December 2022).
  32. Zhang, B.; He, L.; Yi, S.; Ding, H. Multi-resolution analysis of aero-optical effects in a supersonic turbulent boundary layer. Appl. Opt. 2021, 60, 2242–2251. [Google Scholar] [CrossRef]
Figure 1. Schematic diagram of boundary layer density proxy model.
Figure 1. Schematic diagram of boundary layer density proxy model.
Sensors 23 01646 g001
Figure 2. Schematic diagram of simulation threshold of density proxy model.
Figure 2. Schematic diagram of simulation threshold of density proxy model.
Sensors 23 01646 g002
Figure 3. Schematic diagram of ellipsoidal vortex coordinate system.
Figure 3. Schematic diagram of ellipsoidal vortex coordinate system.
Sensors 23 01646 g003
Figure 4. SR of density proxy model changing with control parameters.
Figure 4. SR of density proxy model changing with control parameters.
Sensors 23 01646 g004
Figure 5. The solution results of the hyper-parameters based on Bayesian optimization: (a) iterative scatter distribution of N1, N2; (b) iterative scatter distribution of N3, N4; (c) iteration trajectory of the trust coefficient γ R ; (d) the change curve of loss function.
Figure 5. The solution results of the hyper-parameters based on Bayesian optimization: (a) iterative scatter distribution of N1, N2; (b) iterative scatter distribution of N3, N4; (c) iteration trajectory of the trust coefficient γ R ; (d) the change curve of loss function.
Sensors 23 01646 g005
Figure 6. (a) S R s of the turbulence density proxy model; (b) S R p of the distortion prediction model.
Figure 6. (a) S R s of the turbulence density proxy model; (b) S R p of the distortion prediction model.
Sensors 23 01646 g006
Figure 7. Optical characteristics of XY plane in density proxy model: (a) t = 0   ms ; (b) t = 5   ms ; (c) t = 10   ms ; (d) t = 15   ms .
Figure 7. Optical characteristics of XY plane in density proxy model: (a) t = 0   ms ; (b) t = 5   ms ; (c) t = 10   ms ; (d) t = 15   ms .
Sensors 23 01646 g007
Figure 8. (a) OPD of the turbulence density proxy model; (b) OPD of the flow field by the CFD.
Figure 8. (a) OPD of the turbulence density proxy model; (b) OPD of the flow field by the CFD.
Sensors 23 01646 g008
Figure 9. (a) OPD spectrum distribution of flow direction; (b) OPD spectrum distribution error of flow direction.
Figure 9. (a) OPD spectrum distribution of flow direction; (b) OPD spectrum distribution error of flow direction.
Sensors 23 01646 g009
Figure 10. Monte Carlo verification results under different conditions: (a) distribution of sampling points at 10 km; (b) relative error of OPD at 10 km; (c) distribution of sampling points at 15 km; (d) relative error of OPD at 15 km.
Figure 10. Monte Carlo verification results under different conditions: (a) distribution of sampling points at 10 km; (b) relative error of OPD at 10 km; (c) distribution of sampling points at 15 km; (d) relative error of OPD at 15 km.
Sensors 23 01646 g010
Figure 11. Simulation of distorted star maps: (a) star A at 15 km; (b) star B at 15 km; (c) star A at 5 km; (d) star B at 5 km.
Figure 11. Simulation of distorted star maps: (a) star A at 15 km; (b) star B at 15 km; (c) star A at 5 km; (d) star B at 5 km.
Sensors 23 01646 g011
Figure 12. The simulation results of distorted star maps with the previous study [27].
Figure 12. The simulation results of distorted star maps with the previous study [27].
Sensors 23 01646 g012
Table 1. The gray value at the center of the distorted star maps.
Table 1. The gray value at the center of the distorted star maps.
H / km Center Gray Value of Star ACenter Gray Value of Star B S R p S R s
5143830.06210.0645
102221270.86130.8650
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Yang, B.; Yu, H.; Liu, C.; Wei, X.; Fan, Z.; Miao, J. A Quick Simulation Method for Aero-Optical Effects Based on a Density Proxy Model. Sensors 2023, 23, 1646. https://doi.org/10.3390/s23031646

AMA Style

Yang B, Yu H, Liu C, Wei X, Fan Z, Miao J. A Quick Simulation Method for Aero-Optical Effects Based on a Density Proxy Model. Sensors. 2023; 23(3):1646. https://doi.org/10.3390/s23031646

Chicago/Turabian Style

Yang, Bo, He Yu, Chaofan Liu, Xiang Wei, Zichen Fan, and Jun Miao. 2023. "A Quick Simulation Method for Aero-Optical Effects Based on a Density Proxy Model" Sensors 23, no. 3: 1646. https://doi.org/10.3390/s23031646

APA Style

Yang, B., Yu, H., Liu, C., Wei, X., Fan, Z., & Miao, J. (2023). A Quick Simulation Method for Aero-Optical Effects Based on a Density Proxy Model. Sensors, 23(3), 1646. https://doi.org/10.3390/s23031646

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop