A method based on 3D ray tracing for aero optical wavefront analysis

Optik 126 (2015) 4392–4396

Contents lists available at ScienceDirect

Optik journal homepage: www.elsevier.de/ijleo

A method based on 3D ray tracing for aero optical wavefront analysis Xiaofei Chang a , Tao Wang b,∗ , Shizheng Wan a , Jie Yan a , Wenxing Fu a a b

College of Astronautics, Northwestern Polytechnical University, Xi’an 710072, China R&D Center, China Academy of Launch Vehicle Technology, Beijing 100076, China

a r t i c l e

i n f o

Article history: Received 27 August 2014 Accepted 25 August 2015 Keywords: Aero optics Transmission effect Computational fluid dynamics grids 3D ray tracing Wavefront phase

a b s t r a c t Considering missiles mounted with infrared detection system flying at hypersonic speed, infrared radiation (IR) from the target will suffer light deflection and wavefront distortions when passing through the flowfields surrounding irdome, namely aero-optical effect. The accuracy of detection and identification probability for missile’s IR detectors will be obviously cut down. This paper presents a method based on 3D ray tracing for studying the aero optical light transmission effect existing on the surface of hemispheric irdome. Density distribution surrounding irdome is calculated utilizing Reynolds Averaged Navier–Stokes method and k-epsilon turbulence model. Gladstone–Dale formula is used to establish the relationship between density and refractive index. Generalized regression neural network and Barron operator are used for interpolating refractive index and its gradient respectively so as to obtain the refractive index distribution. Three-dimensional (3D) ray tracing algorithm is derived based on Snell’s Law, which can effectively calculate the optical path difference, Strehl ratio and other parameters under different flight altitudes, speeds and angle of attacks. © 2015 Elsevier GmbH. All rights reserved.

1. Introduction When aircrafts mounted with infrared detection system is high-speed flying in the atmosphere, intense interactions take place between irdome and coming flow, leading to shock wave, expansion wave, turbulent boundary layer and other flow structures. Infrared radiation is strongly affected by air density gradient changes and aerodynamic heating effects in mixing layer when passing through the complex flow field. Wavefront aberration and light deflection arise, resulting in degradation of target image for infrared detection system. All the phenomena are called aerooptical effects [1], which will reduce signal-to-noise ratio of detection system, shorten the target detection range, and descend the probability of target detection and recognition. In one word, aero-optical effects seriously weaken the ability of infrared detection system, making it difficult to detect, identify and track targets. Scholars around the world have kept the research on mechanism, computation theory and measurement technique of aero-optical effects since 1952, when Liepmann [2] first studied turbulent aero-optical effect. In particular, Sutton [3] and Jumper [4], as the representative of many other scholars, studied the basic theory and experimental measurement of optical optics and published a number of documents, leading the direction of the aero-optical

∗ Corresponding author. E-mail address: [email protected] (T. Wang). http://dx.doi.org/10.1016/j.ijleo.2015.08.171 0030-4026/© 2015 Elsevier GmbH. All rights reserved.

research. The study of aero-optical mechanism can be divided into three parts: optical transmission effect, aerothermodynamic effect and heat radiation effect. Among them, optical transmission effect mainly concerns blur, jitter and offset of target image introduced by the dramatic changes of external flow field. Sutton [5] summarized the general equations of aero-optical analysis by researching on the mathematics of average flow field, turbulent boundary layer, free shear layer and non-uniform turbulence. Liu [6] divided aerodynamic flow into two categories: mean flow and turbulent flow, so as to analyze the effect of complex flow field to image degradation. Tromeur [7] used large eddy simulations for studying aero-optical effects. Harri [8] utilized line and surface array wavefront sensor to acquire wavefront aberration directly. Tian [9] concluded the interaction principle between meticulous flow field of supersonic irdome and aero-optical effect through large amounts of wind tunnel tests. Yao [10] also applied Hartmann sensor and shearing interferometer to the measure of optical transmission effect. Higher cost, adjustment inconvenience and measurement errors obviously limit the use of wind tunnel experiment. With the development of computational fluid dynamics (CFD) technology, CFD method is widely used for the study of light transmission and optical imaging through high-speed flowfields, and has become an important research topic owning to its convenience and low cost. Trolinger [11] developed AOTS method. It uses ZEMAX, a well-known ray tracing software, together with experimental data and CFD process to build multi-phase screen function, which can

X. Chang et al. / Optik 126 (2015) 4392–4396

simulate and test aero-optical effect. Wang [12] proposed a modeling strategy for cell transmission of refractive index in hypersonic flowfields based on the CFD grids. In this paper, the parameters of flowfields surrounding hemispheric irdoom under typical flight conditions are obtained through RANS method. Optical transmission effect is analyzed based on 3D ray tracing method. Wavefront phase difference and Strehl ratio are calculated with regard to different scenarios of height, velocity and attack of angle. This subject helps to lay the foundation for hypersonic flowfield light transmission research.

Navier–Stokes equations are utilized to compute flowfields around hemispheric irdome. The differential form of Navier–Stokes equations for compressible viscous fluid in three-dimensional Cartesian coordinate is given by: ∂U ∂FU ∂GU ∂HU ∂Fv ∂Gv ∂Hv + + + + + + =0 ∂t ∂x ∂y ∂z ∂x ∂y ∂z

(1)

where FU , GU , HU are fluxes of inviscid convective vector, and Fv , Gv , Hv are viscid vector fluxes. Four thermodynamic relations are introduced to promote equations p = RT p 1 e= + ui uj 2 ( − 1) = =

C1 T 3/2 T + C2

(2)

Cp R = Pr ( − 1)Pr

where p, e,  and  stands for pressure, total energy per unit mass, kinetic viscosity coefficient and thermal conductivity respectively, R is the gas constant,  is gas specific heat, coefficient C1 = 1.458 × 10−6 kg m−1 s−1 K−1/2 , and Pr is the Prandtl number and usually Pr = 0.72 for perfect gas. Standard k–ε model is used for numerical simulation of supersonic turbulent flow in this paper, where k is for the turbulent kinetic energy and ε is for the turbulent dissipation rate.

n − 1 = KGD 

(4)

For infrared light, the G–D constant is an approximate function of wavelength, and can be expressed by



Based on the density distribution of flowfields near missile irdome, the light propagation principle in the variable refractive index medium, light path related with the incident height and convergence angle of each surface is computed, which is called ray tracing. Therefore, light propagation can be accomplished so that deflection angle, optical path difference, wavefront aberrations, Strehl ratio and other optical information can be calculated. 3.1. Density to refractive index Fermat’s principle describing the path of light transmission shows that the path of light propagation in the flow field depends on the refractive index. Therefore, the refractive index along the transmission path needs to be calculated using the already known flowfield data. Lorentz–Lorenz formula [13] gives the relationship between flowfield density and refractive index as follows: (3)

where KGD is Gladstone–Dale constant,  stands for atmospheric density and n is the refractive index. Even though the atmospheric

1+

7.52 × 10−3 2

 (m3 /kg)

(5)

where  is the wavelength of incident light, and its unit is microns.

3.2. GRNN-based data processing The density of grid points can be converted into corresponding refractive index according to Gladstone–Dale formula given hereinbefore. These grid points are distributed random and nonuniform in space. To get refractive index and its gradient at any point, immediate interpolation with surrounding discrete points is required. What really makes difference is the way of interpolation. Tradition approach is to modify flowfield grids [14]. More specifically, it normalizes grid points to vertices of cuboid at the cost of certain precision. Then, distance-weighted interpolation method is applied to obtain the wanted refractive index. This method is not only process consumption, but also introduces large errors. Proposed by Spechtin 1991, generalized regression neural network (GRNN) is a variation of radial basis function (RBF) network, and also a feed forward neural network model based on the theory of non-linear regression. Benefiting from faster response speed, less training data, GRNN method can easily handle multi-input singleoutput fitting problems. In this paper GRNN approach is chosen to solve the interpolation problem of refractive index field. GRNN uses normal distribution as its probability density function and fitting result is obtained by the following equation.

  n Y exp −Di2 /2 2 i=1 i  2  Y (X) = n 2 i=1

3. CFD-based ray tracing

2 n2 − 1 1 = KGD 3 n2 + 1 

refractive index is highly determined by Mach and Reynolds number, this paper only considers the contribution of atmospheric density to refractive index. Since the refractive index n almost equals to 1, Eq. (3) can be simplified as follows:

KGD ≈ 2.24 × 10−4

2. Flowfields computation

4393

exp −Di /2

(6)

where Di2 = (X − Xi )T · (X − Xi ) is the square of Euclidean distance  from fitting  point x to the ith training point xi . Moreover, exp −Di2 /2 2 is the activate function, representing the contribution weight of ith training point to fitting result. (1) Refractive index interpolation The refractive index of propagation point is obtained as follows. Construct a cube with the unknown point as the center. This cube contains a number of discrete data points and each point corresponding to a refractive index value. Accomplish training of GRNN with spatial coordinate as the input and corresponding refractive index as the output. Expand the cube if the training data is insufficient. Then, the refractive index can be acquired by substituting the coordinate of propagation point to the trained GRNN. (2) Refractive index gradient interpolation Upon completion of interpolating refractive index, the refractive index gradient in three directions also needs to be calculated to complete the ray tracing process. Refractive index gradient for specific mesh node can be solved utilizing the refractive index value of surrounding nodes. Benefiting from its high precision, Barron operator is applied in this article. After interpolating discrete refractive index along XYZ direction, the refractive index gradient of three directions is acquired respectively.

4394

X. Chang et al. / Optik 126 (2015) 4392–4396

Gradient interpolation for point P(i, j, k) is given by

α 4, β 4, γ

1 ∂n(i, j, k) 8 = n(i − 2, j, k) − n(i − 1, j, k) 12 12 ∂x

θ

1 8 ∂n(i, j, k) = n(i, j − 2, k) − n(i, j − 1, k) 12 12 ∂y

1 8 n(i, j, k + 1) − n(i, j, k + 2) 12 12

According to the optical theory, light transmission in the same uniform refractive index follows rectilinear propagation law. However, when it comes to a surface with different refractive indices on both sides, catadioptric phenomenon occurs following law of reflection and law of refraction [15]. Ray in variable refractive index medium will be propagated through multiple refractions only if the spatial position of refracted ray is primarily acquired. A ray in space can be described by a vector with Cartesian coordinate and direction cosine. According to the law of refraction, once known the vector of incident ray and normal line, together with the refractive index on both sides of interface, the refracted ray can be determined. Fig. 1 shows the light transmission through the interface with different refractive indices. (˛1 , ˇ1 ,  1 ), (˛2 , ˇ2 ,  2 ), (˛3 , ˇ3 ,  3 ) and (˛4 , ˇ4 ,  4 ) donate the direction angles for incident ray, refracted ray, normal line and reflected ray respectively. 1 and 2 represent angles of incidence and refraction. n1 and n2 stand for the refractive index on each side of interface. According to the angles formed by incident ray and refracted ray, applying law of refraction, 3D Snell ray tracing can be summarized below

α 3, β 3, γ

3

cos ˛4 = cos ˛1 − 2 cos 1 cos ˛3 cos ˇ4 = cos ˇ1 − 2 cos 1 cos ˇ3

(10)

cos 4 = cos 1 − 2 cos 1 cos 3

3.4. Wavefront analysis From the perspective of modern optics, wave optics is wavefront optics. There is a range of equiphase wave surfaces in wave field. Light field, which reaches the receiving plane, determines the effect of optical transmission. This light field is known as wavefront and can well reflect the type and characteristic of wave. Wavefront properties often derive from optical path length (OPL), which mathematically integrates refractive index along the path of light propagation:



S2

OPL(x, y, t) =

cos ˛2 cos ˇ3 − cos ˇ2 cos ˛3

(11)

n(x, y, z, t)ds S1

sin 2 (cos ˛1 cos ˇ3 − cos ˇ1 cos ˛3 ) sin 1

cos ˇ2 cos 3 − cos 2 cos ˇ3 (8)

sin 2 (cos ˇ1 cos 3 − cos 1 cos ˇ3 ) sin 1

where S1 indicates the starting point of the optical transmission path, S2 indicates the end point, n is the refractive index along the transmission path. Compared to the absolute optical path length, more practical is to obtain optical path difference (OPD): OPD(x, y, t) = OPL(x, y, t) − OPL(x, y, t)

cos 2 cos ˛3 − cos ˛2 cos 3 =

Normal Li ne

Regarding refracted ray as the next incident ray, light propagation is realized step by step. Thus ray trajectory in variable refractive index medium can be calculated. However, critical problem would occur in this procedure if total internal reflection took place. In real field, root cannot be extracted for a negative value, when (n1 /n2 )sin 1 > 1 and the refracted ray becomes non-existent. In this case, the next incident ray should take the reflected ray and its direction cosine is formulated by

cos˛3 cos ˛1 + cos ˇ3 cos ˇ1 + cos 3 cos 1 = cos( 1 − 2 )

=

O

1

2

2

Fig. 1. The principle of transmission through the interface between two different mediums

3.3. 3D Snell ray tracing

=

α 2, β 2, γ

θ

Incident Ray n1 n2 α 1, β 1, γ 1 Interface

(7)

1 8 n(i, j + 1, k) − n(i, j + 2, k) 12 12

1 8 ∂n(i, j, k) = n(i, j, k − 2) − n(i, j, k − 1) 12 12 ∂z +

Refracted Ray

Reflected Ray

8 1 + n(i + 1, j, k) − n(i + 2, j, k) 12 12

+

4

(12)

where OPL(x, y, t) is the mean value of spatial optical path length. The relationship between wavefront phase difference and OPD is computed by the following equation.

sin 2 (cos 1 cos ˛3 − cos ˛1 cos 3 ) sin 1

n1 sin 1 = n2 sin 2

(r, t) = kOPD

(13)

Now the direction cosine of refracted ray is expressed by:



cos ˛2 =

 cos ˇ2 =

 cos 2 =

n1 n2 n1 n2 n1 n2





(cos ˛1 − cos 1 cos ˛3 ) ± cos ˛3



 (cos ˇ1 − cos 1 cos ˇ3 ) ± cos ˇ3



 (cos 1 − cos 1 cos 3 ) ± cos 3





2





2

where k stands for wave number, and k = 2/. The pulse wave of density brings in a statistic phase variance for wavefront, which is always derived from the root mean square value of OPD.

1 − n21 /n22 sin 1 1 − n21 /n22 sin 1



 2

(9)

2

1 − n21 /n2 sin 1

The second term on the right should be positive if 0 < 1 < /2 while negative if /2 < 1 < .

2

ı2 = (kOPDRMS ) =

 2OPD

RMS





2 (14)

where OPDRMS (x, y) = OPD(x, y, t)2 . Strehl ratio (SR) is extensively chosen as the evaluation criteria for an optical system. According to large-aperture approximation,

X. Chang et al. / Optik 126 (2015) 4392–4396

4395

Fig. 2. Wavefront calculation results under different flying conditions.

an estimate of SR using statistic phase variance for imaging system is given by [1]:



SR = exp −ı2



(15)

4. Result analysis Based on the parameters of aerodynamic flow field under different flight conditions, namely height, velocity and angle of attack (AOA), ray tracing is completed to obtain the optical path length. Strehl ratio is computed applying Eqs. (12)–(15) to study the influence of flight conditions to infrared imaging system. Fig. 2 illustrates the OPD distributions under four flight conditions, respectively the height is 6 km, and the velocity is 3.4 Mach number (Ma), and the angle of attack is 0; the height is 10 km, the velocity is 3.4 Ma, and the angle of attack is 0; the height is 6 km, the velocity is 2.6 Ma, and the angle of attack is 0; the height is 6 km, the velocity is 3.4 Ma, and the angle of attack is 5◦ . Strehl ratios for missile borne optical system under different flight conditions are shown in Table 1. From Fig. 2 and Table 1, some conclusions can be drawn: (1) with regard to the same Mach number and angle of attack, the OPDRMS introduced by high-speed flow field decreases rapidly with increasing height. Because the compression effect of density field decreases along with the flying altitude, it leads to smaller changes in OPD and less effect for aero-optical transmission. (2) With regard to the same height and angle of attack, the OPDRMS introduced by high-speed flow field will be increased with the Mach number. The higher Mach number is, the stronger flow field compresses. More dramatic changes in refractive index lead Table 1 Strehl ratio under different flight conditions. Height (km) Mach AOA (◦ ) OPDRMS (m) Strehl ratio

6 3.4 0 1.5548e−7 0.9421

10 3.4 0 0.9668e−7 0.9772

6 2.6 0 1.3953e−7 0.9531

6 3.4 5 1.7390e−7 0.9281

to more optical distortion, which strongly affects the ability of detection system. (3) With regard to the same height and Mach number, the spatial distribution of OPD shifts with angle of attack changing. When the angle of attack is 0, the OPD symmetrically will be distributed along the Z axis. When there is a positive angle of attack, the symmetry axis starts to move along Z axis. From the result of Strehl ratio, attenuation distortion by aero optics becomes more critical when angle of attack exists. 5. Conclusion When missiles mounted with infrared detection system are flying at supersonic velocity, the flowfield surrounding irdome will change dramatically with flight scenarios. Infrared radiation from target will suffer great distortions when passing through the flowfield. As a consequence, offset, blur and other effects are brought into the target images for IR detection systems. In this paper, RANS method and k–ε turbulent model are used to compute the density distributions of flowfields under typical flight conditions. GRNN method is applied for flowfield data interpolation. 3D Snell approach is chosen to fulfill spatial ray tracing. The effect of flight conditions on aero-optical wavefront transmission is analyzed after calculating OPD and Strehl ratio. This paper proposes an efficient way to obtain accurate flowfield data and to calculate OPD that can directly reflect the light transmission effect in complex flowfields. This work not only lays foundation for impact analysis and recovery research on aero-optical degraded image, but also possesses a certain reference value to reveal and correct aero optics. One step further is to analyze the light transmission in time-varying flowfields and establish a high fidelity model for aero-optical light transmission effects. References [1] E.J. Jumper, E.J. Fitzgerald, Recent advances in aero-optics, Prog. Aerosp. Sci. 37 (3) (2001) 299–339. [2] H.W. Liepmann, Deflection and Diffusion of a Light Ray passing Through a Boundary Layer. Report No. SM-14397, Douglas Aircraft Company, 1952.

4396

X. Chang et al. / Optik 126 (2015) 4392–4396

[3] G.W. Sutton, J.E. Pond, R. Snow, Hypersonic Interceptor Performance Evaluation Center: Aero-optics Performance Predictions. AIAA 1993-2675. [4] E.J. Jumper, Fluid-Optic Interactions II (Final report), 22, University of Notre Dame, Notre Dame, Indiana, USA, 2000. [5] G.W. Sutton, E. John, Optics equations for aero-optical analysis, in: International Conference on Applications of Optics and Photonics, Proc. SPIE, 2011, 80012A80012A-8. [6] Chun-sheng Liu, Recent researches on the optical transmission distortions throughout the average flow field, J. Syst. Eng. Electron. 25 (8) (2003) 941–945. [7] E. Tromeur, E. Garnier, S. Pagaut, et al., Large eddy simulations of aero-optical effects in a turbulent boundary layer, J. Turbul. 4 (1) (2003) 5. [8] S.R. Harris, J. Widiker, B. Duncan, Simultaneous MHz rate flow visualization and wavefront sensing for aero-optic, in: 41st AIAA Aerospace Sciences Meeting and Exhibit. January 6–9, Reno, NV, 2003. [9] Lifeng Tian, Mechanism Study of Fine Structures and Aero-optical Effects of Supersonic Flow Around an Optical Dome, National University of Defense Technology, 2011, pp. 4.

[10] Xianghong Yao, Yungang Wu, Yong Chen, Weiming Xie, Experimental study on the optical propagation effect of the flow field surrounding the optical head cover in supersonic wind tunnel, J. Exp. Fluid Mech. 27 (4) (2013) 97–101. [11] J.D. Trolinger, W.C. Rose, Technique for Simulating and Evaluating Aero-optical Effects in Optical Systems. AIAA 2004-471. [12] Tao Wang, Yan Zhao, Dong Xu, Qiuying Yang, Numerical study of evaluating the optical quality of supersonic flow fields, Appl. Opt. 46 (2007) 5545–5551. [13] M. Born, E. Wolf, Principles of optics: electromagnetic theory of propagation, interference and diffraction of light, CUP Arch. (1999) 92–102. [14] R.C. Aguirre, H.J. Catrakis, Aero-optical wavefronts and scale-local characterization in large-Reynolds-number compressible turbulence, AIAA J. 42 (10) (2004) 1982–1990. [15] C. Gomez-Reino, M.V. Perez, C. Bao, Gradient-index Optics: Fundamentals and Applications, Springer Science & Business Media, 2002.