Computation and analysis of backward ray-tracing in aero-optics flow fields

Liang Xu; Deting Xue; Xiaoyi Lv

doi:10.1364/OE.26.000567

1. Introduction

When the vehicle is flying at high speed in the atmosphere, the optical dome of the imaging detection system interacts with the external airflow to form a complex flow field structure, including a shock wave, expansion wave and strong turbulent boundary layer [1,2]. In an active system, the flow field weakens the projected beam energy and deviates it from the target to be irradiated. In the passive system, the image received by the optical system is shifted, blurred and jittered [3]. The influences or effects of this gas flow field on optical propagation and optical imaging are called aero-optical effects [2, 4].

In the past study of aero-optical imaging, most researchers have used forward tracing. The study of aero-optics began in the United States in 80s. The US researchers led by Sutton [5, 6], Clark [7], Jumper [8] and others gradually applied aero-optics research results to the actual vehicle. The THAAD(Terminal High Altitude Area Defense) interceptor, with infrared imaging terminal guidance, and a target-receiving optical side window, is an typical example of aero-optics application. In 1987, Clark and other researchers published an aero-optics performance prediction paper [7]. In this paper, a conical head vehicle like THAAD was modeled. They set up a computation grid and used forward ray-tracing. In 1994, Sutton also published a paper on aero-optics performance prediction for a conical head side window interceptor [6]. In this paper, Sutton also gave a forward tracing scheme, which performed the optical computations along the direction from a defined position in the freestream down to the window. In subsequent studies, scholars basically adopted the forward tracing scheme.

The starting point is required in the forward tracing. Common practice is to pre-define a computational grid before tracing. As shown in Fig. 1, the edge of the grid is selected as the starting line for tracing. The tracing starts from the upper edge of the grid, and ends at the lower edge. The calculation of the entire ray-tracing needs to be done in the preset grid, so the aero-optics flow field needs to be known in advance. The range of the flow field will vary with the height, Mach number, and angle of attack, as shown in Fig. 1. When the incident angle is greater than $θ$ , the grid will not cover some part of the aero-optics flow field, which will cause a tracing failure.

Fig. 1 Illustration of air density distribution and geometry for backward ray-tracing.

Download Full Size | PDF

In this paper the backward ray-tracing scheme is proposed. Performing the traditional forward tracing scheme requires a review of the distribution of the aero-optical flow field beforehand, and then requires preset of a sufficient dense grid to give a sufficiently accurate ray-tracing computation. However, the backward scheme puts a 180° reversion to the tracing procedure and removes the requirement for a preset grid. Then the iterative computation of ray equation in the CFD(Computational Fluid Dynamics) grid is conducted. This approach greatly simplifies the ray-tracing computation. However, a stopping criterion for the backward tracing in the aero-optical flow field is needed to be designed. With this criterion, the disadvantage of paying additional computation to the freestream region is overcome. We observe that the scheme and criterion should bring great benefits: 1) The starting point of tracing is fixed at the center of the entrance pupal of the internal optics. For incident rays of different line-of-sight angles, the starting point is essentially a fixed position, which is no longer a moving position on the edge of the perset grid as in forward tracing; 2) It is not necessary to examine the changes of the aero-optics flow field as is done in forward tracing, then preset a series of starting points that cover all kinds of flow fields for different conditions (And these starting points are usually far from the vehicle); 3) Using a stop tracing criterion, the tracing will end at the outer edge of the aero-optics flow field. This scheme will largely reduce the redundant computations in traditional forward tracing.

2. Calculation of the aero-optics flow field

2.1 Vehicle grid model

In this paper, a typical blunt nosed vehicle is selected as an example. During the flight of the vehicle, the complex flow field structures such as the shock, the expansion wave and the strong turbulence boundary layer vary greatly [9–11]. Therefore, the head mesh quality affects the flow field computations. The closer the flight head is, the denser the grid must be [12–14]. A sufficient grid density relative to an object size is shown in Fig. 2, we set the first level of the boundary of the vehicle head position to a typical value of 0.5 mm, the grid (growth factor) to 1.1, and the whole boundary layer to 4 rows. The mesh was generated in the GAMBIT software [15] and the total grid cells are 295800.

Fig. 2 Local CFD grid.

Download Full Size | PDF

2.2 Flow field calculation

The vehicle flight states chosen are shown in Table 1. It shows a total of $5 \times 6 \times 6 \times 8=1440$ sets of flight states, that we calculated. All the parameters were selected in their typical ranges.

Table 1. Coverage of combinations parameters for which ray deviations were computed.

View Table

The mean flow density contours for three examples are shown in Fig. 3(a)-3(c). The computed Mach number is 3, angle of attack is 0°, and the flight heights are 25km, 15km, 5km, respectively.

Fig. 3 Mean flow density contours.

Download Full Size | PDF

2.3 Gladstone-Dale relation

Gladstone-Dale relation connects the flow field density with the refractive index to obtain a discrete refractive index field, Eq. (1) [12],

n = 1 + K_{G D} ρ,

where, n is refractive index,

ρ

is the density, and

K_{G D}

represents the Gladstone-Dale coefficients, that weakly depend on wavelength. With density expressed in kg/m³ and wavelength expressed in μm,

K_{G D} (λ)

can be approximated by [2]

K_{G D} = 2.23 \times 10^{- 4} (1 + \frac{7.52 \times 10^{- 3}}{λ^{- 2}}),

where

λ

is the light wavelength in μm, and

K_{G D}

is in m³/kg.

2.4 Quadrilateral mesh interpolation

Assuming that the CFD grid cells, on the x-y physical plane, are as shown in Fig. 4. The four vertices of the quadrilateral are 1, 2, 3 and 4. The variables u₁, u₂, u₃ and u₄ represent refractive index. Then the refractive index value u, which falls in this quadrilateral, can be expressed by bilinear interpolation as [15, 16]:

u = u_{1} φ_{1} + u_{2} φ_{2} + u_{3} φ_{3} + u_{4} φ_{4}

Where

φ_{i}

(

i

= 1, 2, 3, 4) can be calculated by the following formula:

{\begin{cases} φ_{1} = \frac{1}{4} (1 - ξ) (1 - η) \\ φ_{2} = \frac{1}{4} (1 + ξ) (1 - η) \\ φ_{3} = \frac{1}{4} (1 + ξ) (1 + η), \\ φ_{4} = \frac{1}{4} (1 - ξ) (1 + η) \end{cases}

where

ξ

and η are the inner coordinates of the quadrilateral, and there are also the computational coordinates for bilinear interpolation shown in Fig. 4(b).

Fig. 4 Refractive index bilinear interpolation.

Download Full Size | PDF

The relationship between the physical coordinate system and the computational coordinate system of bilinear interpolation is given by Eqs. (5):

{\begin{cases} x = \sum_{i = 1}^{4} φ_{i} (ξ, η) x_{i} \\ y = \sum_{i = 1}^{4} φ_{i} (ξ, η) y_{i} \end{cases}

where

x_{i}

and

y_{i}

are the coordinate values of the i vertices of the quadrilateral in the physical plane.

Substituting Eqs. (4) into Eqs. (5), we obtain the functions of ξ and η about x and y, then putting the coordinate values of the tracing point $(x, y)$ into the obtained functions, we get the coordinate values in the ξ-η plane. Then we will substitute the coordinate values of the ξ-η plane into Eqs. (4) to get the values $φ_{i}$ . Finally, we put the $φ_{i}$ into the Eq. (3) and get the refractive index u.

3. Backward ray-tracing computation

3.1 Backward tracing

Figure 5 illustrates the ray-tracing process from the (k-1)th step to the kth step. Assuming the contour line equation at the (k-1)th step is: $y_{k + 1} = f_{k - 1} (x)$ , at the point $(x_{k - 1}, y_{k - 1})$ , the tangent equation is

y - y_{k - 1} = f_{k - 1}^{'} (x_{k - 1}) (x - x_{k - 1})

The angle between the tangential line and the horizontal line is

Fig. 5 Ray-tracing process from the (k-1)th step to the kth step.

Download Full Size | PDF

γ_{k - 1} = {_{π + arc \tan f_{k - 1}^{'} (x_{k - 1}), a r c \tan f_{k - 1}^{'} (x_{k - 1}) < 0}^{arc \tan t f_{k - 1}^{'} (x_{k - 1}), arc \tan f_{k - 1}^{'} (x_{k - 1}) \geq 0}

Generally consider $f_{k - 1}^{'} (x_{k - 1}) \neq 0$ , the equation for the normal is

y - y_{k - 1} = \frac{1}{f_{k - 1}^{'} (x_{k - 1})} (x - x_{k - 1})

Set the (k-1)th step of the incident ray expression

l_{k - 1}

as:

a_{k - 1} x + b_{k - 1} y + c_{k - 1} = 0

,

a_{k - 1} \neq 0

, Then the incident angle is

θ_{k - 1} = arc \cos \frac{(1, f_{k - 1}^{'} (x_{k - 1})) \cdot (a_{k - 1}, b_{k - 1})}{| 1, f_{k - 1}^{'} (x_{k - 1}) | | a_{k - 1}, b_{k - 1} |}

Suppose

β_{k - 1}

is the emergence angle, then Snell's Law gives

\frac{\sin θ_{k - 1}}{\sin β_{k - 1}} = \frac{n_{k}}{n_{k - 1}}

It gives the angle between the exit ray and the horizontal plane as

β_{k - 1} + γ_{k - 1} - \frac{π}{2}

at the point

(x_{k - 1}, y_{k - 1})

. Then the expression for the kth step of the incident light is

y - y_{k - 1} = \tan [β_{k - 1} + γ_{k - 1} - (\frac{π}{2})] (x - x_{k - 1})

Combining Eq. (11) with the contour line of the kth step $y_{k} = f_{k} (x)$ givens the intersection $(x_{k}, y_{k})$ of the two lines. So the light transmission path is from $(x_{k - 1}, y_{k - 1})$ to $(x_{k}, y_{k})$ [17].

3.2 Backward ray-tracing stopping criterion

The end points of the ray-tracing can be determined by stopping criterion that allows ray-tracing stop at the outer boundary of the aero-optics flow field. The criteria for stopping guidelines can be defined by

| n_{p = q} - n_{f} | \leq δ, q = 0, 1, \dots, m,

where,

n_{p + q}

is refractive index at the (p + q)th step,

n_{f}

is the refractive index at the local freestream,

δ

is a small positive number. The value of delta is determined by the free flow of refractive index fluctuations, and can usually be taken as 10⁻⁴. The variable

m

represents the number of additional traces, and it is used to ensure that the ray-tracing has indeed reached the boundary of the aero optical flow field. The value of the positive integer

m

is determined by the advancing geometry step

Δ s

of the ray-tracing, and can be taken as

m = [0.003 / Δ s] + 1

. For example, if

Δ s

is 0.001 then

m

is 4, The ray-tracing stopping criterion is that

| n_{p} - n_{f} |

,

| n_{p + 1} - n_{f} |

,

| n_{p + 2} - n_{f} |

,

| n_{p + 3} - n_{f} |

and

| n_{p + 4} - n_{f} |

are all less than

δ

.

4. Experimental arrangement and results analysis

The variables Alt, Ma, AoA and LOS are the height, Mach number, angle of attack and line-of-sight angle, respectively. The ray-tracing point starts at the internal sensor behind the window, and the tracing step is 0.8mm. The line-of-sight angle is defined as the angle between the light and the vehicle’s body axis.

In order to facilitate the analysis, the ray-tracing is divided into three sections in turn: 1)A-B, tracing within the vehicle’s body, 2)B-C, tracing in the non-uniform flow field, which is also the aero-optics flow field, and 3)C-D, tracing in the freestream. In the backward tracing, the tracing direction is from A to D, so the point A is the starting point that locates inside the body, point B is on the surface of the vehicle’s body, point C is on the boundary of the aero-optics flow field, and the point D is the ending point that locates in the freestream. While in the forward tracing, point A is the ending point and point D is the starting point. Note that in each direction, the point A is a fixed point that locates at the center of the entrance pupal of the internal optics.

Because the refractive index fluctuation in the freestream is very little compared with that in the aero-optics flow field, the light propagation path in the freestream is regarded as a straight line within the working radius of the airborne optical systems. So the ray-tracing computations in the freestream are redundant, which means the shorter the C-D section, the lesser the redundant computations.

Figure 6 shows the backward ray-tracing results at LOS angles from 5° to 85°, with a step length of 20°. Taking the 85° LOS angle for example, the points A, B, C and D are marked on the refractive index curve, and they are also projected onto the X-Y plane for clearly. The ray-tracing starts from point A and ends at point D. In Fig. 6, the point D is very close to point C, which means the tracing ends just at the outer edge of the aero-optics flow field. So the redundant computations are very less in the backward tracing. In this study, the parameter m in the stopping criterion is 4, so the tracing computations in the freestream are only 4 steps. Along the 85° LOS angle tracing curve, from the direction of A to D, the refractive index increases sharply at point B, and it decreases rapidly at point C. This is because the air density in the aero-optics flow field is much higher than that in the vehicle’s inner body and in the freestream. Comparing the different LOS angle tracing curves, at point B, the refractive index becomes higher as the LOS angle decreases. However, at point C, the refractive index varies slowly with different LOS angles. These variations are consistent with the mean density contour distributions of the flow fields shown in Fig. 3.

Fig. 6 Backward Ray-Tracing. Example of the special distribution of refractive index in front of the vehicle optics. Plots are for various angles LOS from the vehicle axis. (Alt = 25km, Ma = 3, AoA = 0°)

Download Full Size | PDF

Figure 7 shows another set of backward ray-tracing results for different flying altitudes range from 5 km to 25 km. Investigating the points A-D marked on the 5 km curve can draw the same conclusion as that in Fig. 6. However, the refractive indexes of points A-D of the 15 km and 25 km curves are lower than that of the 5 km case, which are consistent with the density variation.

Fig. 7 Backward Ray-Tracing. Example of the special distribution of refractive index in front of the vehicle optics. Plots are for various altitudes. (Ma = 3, AoA = 0°, LOS = 5°).

Download Full Size | PDF

In the forward scheme, the starting points are not at a fixed position, and they are usually set to be on a regular curve such as a circular or a line in the freestream far from the vehicle. Taking the circular case for example, Fig. 8 shows the forward tracing results under the same conditions arranged for the forward case of Fig. 6. The points A, B, C and D are marked on the 5 km curve as an example. Differing from the backward cases shown in Fig. 6 and 7, the forward tracing starts at the point D and ends at the point A. It is seen that the section C-D in Fig. 8 is much longer than that in Fig. 6, with compared to other sections A-B and B-C under the same scale. It means that much more redundant tracings will be paid for section C-D in the forward tracing than that in the backward tracing. Figure 9 gives the quantitative comparison between the tracing steps of the two schemes.

Fig. 8 Forward Ray Tracing. Example of the special distribution of refractive index in front of the vehicle optics. Plots are for various angles LOS from the vehicle axis. (Alt = 25km, Ma = 3, AoA = 0°).

Download Full Size | PDF

Fig. 9 Comparison of forward tracing steps and backward tracing steps (Alt = 25km, Ma = 3, AoA = 0°, LOS = 5°).

Download Full Size | PDF

Transformation to the actual, forward ray requires a distant, starting point for it, but the non-uniform part of the flow field is within the range of the forward ray-tracing. As shown in Fig. 8, the starting points for forward tracing is set at D, the forward ray-tracing can accurately trace, but we need to pay attention to two points: First, consider the starting point D. When the ray is traced in different LOS angles, it is necessary to change the starting point, so as to reach the fixed ending point A, which locate at the center of the entrance pupal of the internal optics. Second, it is redundant to calculate ray-tracing in the freestream area. In the process of each direction of ray-tracing, the tracing range must include the non-uniform flow fields. But in different flight conditions, the non-uniform flow fields are different. In order to guarantee the coverage of the preset grid for forward tracing to all of the propagation paths of incident rays in the non-uniform flow fields, the grid would have to be sufficient large. So the range of tracing would increase a lot.

In this study, two kinds of starting point trajectory were designed for the forward tracing. Figure 8 has shown the circular trajectory, whose equation is expressed by ${(x - 0.08)}^{2} + y^{2} = {0.3}^{2}$ . In order to give a further illustration, the tracing steps of a straight line trajectory were also shown. The straight line is expressed by $35 x - 27 y + 4.55 = 0$ . Comparing the three curves in Fig. 9, it is seen that the backward tracing has the least tracing steps. At LOS angle 5°, the backward tracing can save about 80% computations compared with the circular starting point trajectory. At LOS angle 45°, the backward tracing steps is almost as much as that of the straight trajectory. This is because the straight line is carefully designed according to the flow field distribution, which has been computed and drawn beforehand. Let the line just covers the aero-optics flow field.

In the backward scheme, the starting position is a fixed point, not a moving point as in the forward scheme, which will also reduce the computational errors. In Figs. 6 and 7, the ray-tracing sections C-D have used the stop-tracing criterion, which make the tracing stop at the outer edge of the aero-optics flow field adaptively. So the redundant calculations payed for the large grid have been avoided.

The correctness of the backward tracing results is guaranteed by the principle of optical path reversibility in geometrical optics. In our project's passive aero-optics imaging study, it is the radiation emission from the target's exhaust or reflection of sunlight from the target's body by a distant target that arrives at the sensor. In the project's active, laser illumination aero-optics study, the light emitted from the laser on our interceptor, and arrives at the target. Therefore, the backward tracing used in imaging study will turn into a forward tracing in our laser illumination study. In such a system, the backward tracing method and the stop-tracing criterion proposed in this paper are also applicable.

5. Conclusions

In this paper we propose a backward ray-tracing method for aero-optics simulations. Differing from previous forward tracing method that performed tracing in a preset grid, we discard the preset grid and reverse the tracing procedure, videlicet, start tracing from the sensor and end it at the free stream boundary. As compared to the forward way, the advantage of this backward method is that it performs tracing in an open direction facing to the free stream, no longer a closed direction facing to the sensor. So the backward method does not need to preset the scope of the tracing nor analyze the non-uniform flow-field distributions beforehand. This so as to guarantees the coverage of a much smaller preset grid to all propagation paths of incident rays in the non-uniform flow field. It allows us to present a criterion that ensures that the tracing stop at the outer bounder of the aero-optical flow field automatically. The backward method combined with stopping criterion greatly simplify the ray-tracing computations in aero-optics simulations, and also avoids the redundant computations that are inherent in forward tracing within the free stream region.

Funding

National Natural Science Foundation of China (NSFC) (61308120, 61463029); Reserve Talents Project of National High-level Personnel of Special Support Program (QN2016YX0324); Reserve Talents Project of National High-level Personnel of Special Support Program (Xinjiang[2014]22).

References and links

1. X. L. Yin, “A new subdiscipline of contemporary optics–aero-optics,” China Engineering Science 7(12), 1–6 (2005).

2. X. L. Yin, Principle of Aero-optics (China Aerospace Publishing House, 2003).

3. L. Xu and Y. Cai, “Influence of altitude on aero-optic imaging deviation,” Appl. Opt. 50(18), 2949–2957 (2011). [PubMed]

4. G. C. Li, Aero-optics (National Defense Industry Press, 2006).

5. G. W. Sutton, J. E. Pond, and R. Snow, “Hypersonic interceptor performance evaluation center aero-optics performance predictions,” presented at the 2nd Annual AIAA SDIO Interceptor Technology Conference, Albuquerque, NM , America (NY)6–9(June), (1993).

6. G. W. Sutton, J. E. Pond, and R. Snow, “Hypersonic interceptor aero-optics performance predictions,” J. Spacecr. Rockets 31(4), 592–599 (1994).

7. R. L. Clark and R. C. Ferris, “A numerical method to predict aero-optical performance in hypersonic flight,” presented at the AIAA 19th Fluid Dynamics Plasma Dynamics and Lasers Conference, Honolulu, Hawaii, 8–10 June (1987)

8. E. J. Jumper and R. J. Hugo, “Optical phase distortion due to turbulent-fluid density fields: Quantification using the small aperture beam technique,” presented at the AIAA 23rd Plasmadynamics & Lasers Conference, Nashville, TN , America (NY)6–8(July), (1992).

9. S. Z. Wan, X. F. Chang, and J. Yan, “Research method of aero-optical transmission effect for IR air to air missiles,” Journal of Northwestern Polytechnical University 33(04), 621–626 (2015).

10. G. Guo, H. Liu, and B. Zhang, “Aero-optical effects of an optical seeker with a supersonic jet for hypersonic vehicles in near space,” Appl. Opt. 55(17), 4741–4751 (2016). [PubMed]

11. T. Wang, Y. Zhao, D. Xu, and Q. Yang, “Numerical study of evaluating the optical quality of supersonic flow fields,” Appl. Opt. 46(23), 5545–5551 (2007). [PubMed]

12. W. Merzkirch, Flow Visualization, 2nd ed,(Academic, 1987).

13. L. Xu and Y. L. Cai, “Imaging deviation through non-uniform flow fields around high-speed flying vehicles,” Optik – Int. J. Light Electron Opt 123, 1177–1182 (2012).

14. T. Wang, Y. Zhao, D. Xu, and Q. Y. Yang, “CFD grids-based transmission model of the rays propagating through the hypersonic flow field,” Acta Armamentarii 29(03), 282–286 (2008).

15. C. Yan, Study on Computational Fluid Dynamics And Its Application (Beihang University, 2006).

16. B. Z. Zhang, J. A. Yin, and H. J. Zhang, Numerical Methods of Fluid Mechanics. (China Machine Press, 2003).

17. W. X. Yang, C. Chao, M. Y. Ding, and C. P. Zhou, “Analysis of Pneumatic Optical Effect of Turbulent Flow Field in Supersonic / Hypersonic Vehicle,” Photosystems 36(01), 88–92 (2009).

Parameter	Ranges	Intervals
Vehicle Altitude	0-25km	5km
Mach number	0.5-3M	0.5M
Angle of attack	0-60°	5°,15°
Line-of-sight angle	5-85°	10°

Computation and analysis of backward ray-tracing in aero-optics flow fields

Abstract

1. Introduction

2. Calculation of the aero-optics flow field

2.1 Vehicle grid model

2.2 Flow field calculation

2.3 Gladstone-Dale relation

2.4 Quadrilateral mesh interpolation

3. Backward ray-tracing computation

3.1 Backward tracing

3.2 Backward ray-tracing stopping criterion

4. Experimental arrangement and results analysis

5. Conclusions

Funding

References and links

Cited By

Figures (9)

Tables (1)

Equations (12)

Optics Express