Measurement Techniques, Vol. 52, No. 5, 2009

ALGORITHM FOR CALCULATING THE REFRACTION PATTERNS OF A PLANAR LASER BEAM IN AN OPTICALLY INHOMOGENEOUS MEDIUM

K. M. Lapitskii, I. L. Raskovskaya, and B. S. Rinkevicius

UDC 535.31:681.7.001

An algorithm for simulation of the distribution of a planar laser beam in a three-dimensional, optically inhomogeneous medium is proposed. A technique for numerical calculation of the trajectories of the beam in a geometrically optical approximation is described. A calculation of refraction patterns in probing of a thermal boundary layer near a cylindrically shaped heated body that takes into account edge effects is performed. Key words: refraction, laser beam, optically inhomogeneous medium, boundary layer, temperature inhomogeneity, edge effects.

The broad application of laser methods of diagnostics of flows of liquids and gases derives from their basic advantages as compared with other methods. Through the use of laser radiation it becomes possible to avoid having to introduce distortions into the particular medium and also assures that the measurements will be noncontact, inertialess, and remote [1]. The principal requirement imposed on the particular medium is that it must be optically transparent. The characteristics of the distribution of the index of refraction of the medium that are obtained in an experiment may subsequently be converted into the distribution of other physical quantities (temperature, density, salinity, etc.). Natural convection in a liquid near heated bodies is one of the more widespread types of flows. The principal difficulty in the study of such flows is related to the three-dimensional nature of the flow, along with their substantial irregularity and the fact that they are non-steady-state in nature. Numerical simulation on the basis of modern mathematical models of convective flows requires significant computational resources [2]. Such flows constitute a classical form of optically inhomogeneous media, since the index of refraction of the liquid depends on temperature. A number of different optical methods [1], in particular, the method of laser refractography [3], have been developed for use in the study of such flows. An analytic solution of both the direct problem of calculating the trajectory of the propagation of geometric-optical rays in an optically inhomogeneous medium as well as the inverse problem of reconstruction of the law of distribution of the index of refraction of the medium is possible only for axisymmetric bodies [4] (for example, spheres, ellipsoids of rotation, cylinders) or bodies of planar form, disregarding edge effects. In all other cases, such calculation can be performed only by means of numerical methods. A computer-based laser refraction method [3] has been developed for purposes of diagnostics of optically inhomogeneous flows of liquids and gases. In this method, a flow is probed by a plane laser beam, i.e., an astigmatic beam with elliptical cross-section the dimension of which is substantially greater along one of the axes than along another axis; the plane has also been called a “laser plane” [5]. Edge effects near the edges of bodies of finite dimension [6] as well as the variation in

Moscow Power Engineering Institute (Technical University), Moscow, Russia; e-mail: [email protected]. Translated from Izmeritel’naya Tekhnika, No. 5, pp. 36–39, May, 2008. Original article submitted February 26, 2009. 494

0543-1972/09/5205-0494 ©2009 Springer Science+Business Media, Inc.

a

b

Fig. 1. Geometric parameters of problem: a) trajectory of beam with unit vectors of tracking trihedron I, v, b in field of gradient of index of refraction; b) rotation dϕ of tangent vector of beam in osculating plane α.

the configuration of a regular temperature field [7] lead to a substantial variation in the trajectory of the propagation of the radiation and, as a consequence, in the form of the laser plane. Through the use of algorithms for numerical calculation of the trajectory of beams in an optically inhomogeneous medium and the newly developed computer program that is presented in the article, it will become possible to obtain refraction patterns (tracks of the laser plane at the exit from an optical inhomogeneity) for cases that correspond to physically performed measurements. 1. Derivation of computational relationships. The trajectory of a beam with unit vectors of the tracking trihedron I, v, b along the ray, where I is the tangential unit vector, v the unit vector of the principal normal, and b the unit vector of the trajectory binormal, is depicted in Fig. 1a. The two vectors I and v lie in what is known as the osculating plane α (Fig. 1b) and, as is shown in [4], are connected by the relationship dI/ds = Kv, where the differential of an arc of the trajectory ds = ( dx )2 + ( dy )2 + ( dz )2 , and the curvature of the beam K = (gradn/ nsinϕ), ϕ is the acute angle between the unit vector I and the vector grad n. The binormal unit vector b may be determined for known I and v as b = [l, v]. The unit vectors v and b lie in the plane β (cf. Fig. 1b) perpendicular to the tangent vector I. If the dependence of the index of refraction on coordinates n(x, y, z) is known, a well-known differential equation [8] may be used to find the equation of the trajectory of the beam. By means of this equation, it becomes possible to determine the unit vector I(x, y, z) tangent to the beam at every point of its trajectory: d[n(x, y, z)l(x, y, z)]/ds = gradn.

(1)

Differentiating the left-hand side of (1), we obtain dl/ds = (gradn)/n – (l/n)(dn/ds).

(2)

495

The quantity on the right-hand side of Eq. (2) constitutes the transverse gradient [4] the direction of which coincides with the principal normal vector v. The purpose of this algorithm is to determine the trajectory of a beam; this may be done directly from (2) if the boundary value of the vector I0 is given and with stepwise numerical integration with respect to ds. In Fig. 1a, values of the unit vectors with subscript 1 correspond to the next step. However, it is practically more convenient to pass from the vector equation (2) to a scalar equation, in which the turning angle dϕ of the tangent vector I relative to the direction gradn on the trajectory segment ds is the desired function (cf. Fig. 1b). For this purpose, we expand the vector I with respect to the two unit vectors: p, which coincides with the direction of the vector grad n, and q, a vector perpendicular to gradn, while the unit vectors I, p, and q all lie in the plane l = pcosϕ + qsinϕ. We represent the left-hand side of (2) in the form dl/ds = (dl/dϕ)(dϕ/ds),

(3)

dl/dϕ = –psinϕ + qcosϕ = v.

(4)

where

Substituting (3) and (4) into (2) and performing scalar multiplication of the left- and right-hand sides of the resulting equation by v, we obtain dϕ/ds = (vgrad n)/n = (grad n/n)sinϕ. (5) From (5) it follows that dϕ = (grad n/ n) sin ϕ ( dx )2 + ( dy )2 + ( dz )2 .

(6)

We determine the vector I1, which is obtained by rotation by the angle dϕ in the plane grad n (see Fig. 1b), as l1 = lcosdϕ + vsindϕ.

(7)

v = [gradn – l(dn/ds)]/(gradnsin ϕ),

(8)

dn/ds = (dn/dr)(dr/ds) = lgrad n =gradncosϕ.

(9)

Hence, by (2),

where

For given directional cosines of the vectors l = l0(cosα0, cosβ0, cosγ0) and grad0 n, as well the value of n0 at the entrance to the optically inhomogeneous medium at the point (x0, y0, z0), the values of ϕ and dϕ and the tangent vector I1 corresponding to the first step along the trajectory arc ds may be found from (6)–(8) (Fig. 1a). The construction of the algorithm that will be described below is based on quantization of the computational region and step-by-step determination of the beam trajectory. 2. Algorithm for numerical calculation of the trajectory of beams. Without loss of generality, we will assume that the direction of propagation of the radiation coincides with the Z-axis. We partition the computational region into cells of dimension ∆xk, ∆yk, and ∆zk within which the value of the index of refraction will be assumed to be constant, i.e., nk = n(xk, yk, zk). Therefore, the segment of the trajectory within each cell constitutes a straight line. The dimensions of the cells of the computational grid may be selected on the basis of a comparison of a given solution with the results of analytic calculations for test problems. 496

As the boundary condition, it is first necessary to specify the tangent unit vector to the beam trajectory I0 by means of its directional cosines (cosα0, cosβ0, cosγ0) and its point of entry into the optically inhomogeneous medium, which has the coordinates (x0, y0, z0). Let us consider an arbitrary kth cell of the computational region. The gradient of the index of refraction is determined by the relationship grad k n = grad n

x = xk , y = yk , z = zk

.

Once Ik is known, we calculate the angle ϕk between the unit vector Ik and the vector gradk n on the basis of the scalar product: cos ϕ k =

grad kx n cos α k + grad ky n cos β k + grad kz n cos γ k grad k n

.

Substituting in (6) finite increments for the differentials, we obtain an expression for the turning angle ∆ϕk of the vector Ik in a cell: ∆ϕ k = (grad k nk/ nk ) sin ϕ k (∆x k )2 + (∆y k )2 + (∆z k )2 .

We calculate the vector vk on the basis of (8) and (9), after which we determine the vector Ik+1 in the succeeding cell from (7) thus: lk+1 = lk cos∆ϕk + vk sin∆ϕk. As a result, we obtain the coordinates (xex, yex, zex) and the directional cosines (cosαex, cosβex, cosγex) of the beam vector Iex. The refraction pattern at the exit from the optically inhomogeneous medium may be calculated by specification of a family of beams at the entry to the inhomogeneity. 3. Testing the algorithm used to calculate the trajectory of a light beam. The RELAS (REfraction of LAser Sheet) program, which may be used to solve a problem related to determining the propagation of a laser plane in a spherically stratified, optically inhomogeneous medium, was developed for the purpose of testing the algorithm that has been considered here. The application was created in a Delphi programming medium and may function on computers with the Windows operating system. The following system requirements are recommended: 64 Mbyte internal memory, 10 Mbyte free disk space. Through the use of the program it becomes possible to calculate refraction patterns in spherically symmetric, optically homogeneous media with temperature distribution according to an exponential or Gaussian law with output of the results in the form of graphics and with output to a test file. The computation time for a single refraction pattern on a Pentium III computer with clock frequency 600 MHz is in the range 1–2 min, depending on the number of grid cells. Let us consider a previously solved problem in the following formulation, which approaches an experiment. A heated cylinder with semispherical base of radius R = 17 mm is placed in a vessel of cubic form with length of an edge D = 100 mm. The vessel is filled with water, while the lateral and lower walls of the vessel are kept at a constant temperature. A problem of numerical simulation of the temperature field near a heated object with natural convection was considered in detail in [9], and from the results of the calculation, the structure of the temperature field beneath the base of this cylinder was found to be close to the spherical symmetric. We will consider the case of propagation of radiation beneath a cylinder, where the temperature field is spherically symmetric with a sufficient degree of precision. We arrange for the coordinate origin to coincide with the center of the hemisphere and specify the temperature distribution for a distance r ≥ R with respect to the radial coordinate of an exponential dependence of the form T(r) = T0 + ∆Texp[–(r – ∆R – R)2/a2],

(10)

where T0, ∆T, ∆R, and a are the parameters of the temperature field model [5]; r is the radial coordinate, r = (x2 + y2 + z2)1/2; and R is the radius of the cylinder. The parameter T0, is determined by the temperature of the walls of the vessel containing 497

Fig. 2. Calculation of refraction patterns for passage beneath a heated cylinder with hemispherical base in water: 1) NZ = 100; 2) NZ = 1000; 3) NZ = 10000 and analytic calculation.

the liquid, and the value of T(r = R) is equal to the temperature at the surface of the heated body while the ratio ∆T/a corresponds to the mean gradient of the temperature field in a boundary layer of thickness a determined from the level of e-fold drop in the temperature. The temperature gradient with r = R is characterized by a shift ∆R. To simplify the discussion, the axis of rotation X is directed downwards. Let us consider the case with a temperature distribution of the form (10) with the following parameters: T0 = 18°C; ∆T = 60°C; ∆R = –0.14 mm; R = 17 mm; a = 0.70 mm.

(11)

The dependence of the index of refraction of water on temprtuer for a wavelength λ = 0.6328 µm is found from the formula n(T) = 1.3328 – 0.000051T – 0.0000011t 2 [10]. Suppose that a laser plane 20 mm in width parallel to the Y-axis and having an initial coordinate x0 = R + 0.1 = 17.1 mm propagates along the Z-axis. In order to conduct the calculation according to the technique described in sec. 1, we represent the laser plane as a set of beams with the following initial conditions: x0i = 17.1 mm; y0i = –10 + 0.5i mm; z0i = –50 mm; α0i = β0i = 0; i = 0, 1, ..., 20. The calculations halt once a value of the coordinate zi = 50 mm is attained. The calculation was performed for three values of the number of partitions NZ of the computational region along the Z-axis of 100, 1000, and 10000, i.e., a step along the Z-axis equal to 1, 0.1, and 0.01 mm, respectively. An analytic calculation of the refraction patterns in the indicated optically inhomogeneous medium according to the technique described in [5] was also performed. The results of the calculation are presented in Fig. 2. Curves 1, 2, and 3 correspond to NZ = 100, 1000, and 10000. It follows from the graphs that the refraction pattern obtained from the numerical algorithm approaches that calculated analytically as the dimensions of the cell of the computation region are reduced. For curve 3, the relative error of the calculation according to the analytic technique does not exceed 0.5%. To achieve this level error, the required minimal dimension of a cell of the computation region must not exceed 10 µm, or a/70 times the thickness of the boundary layer. Note that in this calculation we selected a sufficiently small dimension of the computation region along the X- and Z-axes. In typical experiments, the length of an optical inhomogeneity along the Z-axis reaches 200–300 mm, i.e., the computation error will grow while the dimension of the cell remains unchanged. In this case, it is recommended that an adaptive grid with very small step in a region of large temperature gradients (in the boundary region) and a large step in regions where the temperature gradient is close to zero, be selected. 4. Results of simulation of the trajectory of a light beam in a cylindrically inhomogeneous medium. In practice, the problem of studying a boundary layer near cylindrically shaped bodies placed in a liquid is of critical importance. The most substantial decreases in the thermophysical properties occur near the line connecting the base with the lateral surface. 498

Fig. 3. Calculated 1 and experimental 2 refraction patterns for passage beneath the base of a heated cylinder.

Fig. 4. Computational 1 and experimental 2 refraction patterns with passage near the lateral surface of a heated cylinder in water.

The resulting edge effects introduce a major contribution to the flow irregularity, which makes it difficult to study the bottom region by means of refractometric methods. A solution of the inverse problem, which involves reconstruction of the temperature field from an already obtained refraction pattern, is impossible to solve analytically, though an analysis may be performed only by a comparison of the computational and experimental refraction patterns. Let us consider the problem in the following formulation. A heated cylinder with radius R = 17 mm and height H >> R and with base which, unlike the experiment in sec. 3, is planar, is placed in a vessel of cubic shape with length of an edge D = 100 mm. The vessel is filled with water and the lateral and lower walls of the vessel are kept at a constant temperature. We place the coordinate origin at the center of the base of the cylinder (the X-axis is the axis of rotation and the Z-axis the axis of propagation of the laser plane; the X-, Y-, and Z-axes form a right-handed coordinate system) and specify the temperature distribution to be an exponential dependence of the (10) form with parameters (11) as follows: beneath the base of the cylinder the temperature inhomogeneity is plane-stratified, while near the lateral surface it is spherically stratified. The form of the temperature field in an XZ section corresponds to a problem that was considered in [6]. Let the laser plane have a width of 20 mm and suppose it is parallel to the Y-axis and have an initial coordinate x0 = R + 0.1 = 17.1 mm and propagate along the Z-axis. In order to perform the calculation by the technique described in Sec. 1, we represent the laser plane by a set of beams with the following initial conditions: 1) z0 = –50 mm; x0 = 0.1 mm; y0 varies from –30 to 30 mm with the following adaptive step. For y< 16 mm and y> 18 mm, y0 varies every 2 mm, while in the region 16 <y< 18 mm, every 0.1 mm; α0 = β0 = 0; this case corresponds to passage of the laser plane beneath the base of the cylinder; 2) z0 = –50 mm; y0 = R + 0.1 mm = 17.1 mm; x0 varies from –5 to 9 mm with the following adaptive step. For x < –1 mm and x > 1 mm, x0 varies every 1 mm, while in the region –1 < x < 1 mm, every 0.1 mm; α0 = β0 = 0; this case 499

corresponds to passage of the laser plane parallel to the lateral surface of the cylinder. The calculations were halted once a value of the coordinate z = 50 mm was achieved. The calculation was performed with NZ = 1000, i.e., the step along the Z-axis was equal to 0.1 mm. Computational as well as experimental refraction patterns that coincided with the former are presented in Figs. 3 and 4. An analysis of these patterns demonstrates a basically high-qualtiy coincidence of the results of the calculation with experiment. The small difference along some segments is explained by incomplete correspondence of the constructed model of the edge effects to the experiment, and this should be improved. To arrive at a qualitative estimator, a calculation for known parameters of the temperature field taken from an experiment must be performed. Thus, the newly developed algorithm for calculating the propagation of an astigmatic laser beam in a three-dimensional optically inhomogeneous medium is based on stepwise calculation of individual segments of the trajectory in an approximation to a constant value of the index of refraction within a single cell. A computer program for implementation of the algorithm was created. A calculation of refraction patterns in a spherically stratified temperature field from numerical and analytic algorithms shows that the dimension of the cell must be at least a/70 to achieve a relative error in the results of the calculation of 0.5%. The applicability of the present algorithm for the purpose of calculating the trajectory of beams in order to analyze optically inhomogeneous media with edge effects was shown as an example of the propagation of a laser plane near a heated cylinder placed in a cold fluid. The computed refraction patterns are basically similar to the experimental patterns, and the small difference points to incomplete correspondence of the constructed model of the edge effects to the experiment. More precise validation of the results would entail an improvement in the model of the edge effects based on experimental data and more precise thermophysical calculations. The present study was completed with the support of the Russian Foundation for Basic Research (Grant No. 7-07135229).

REFERENCES 1. 2. 3. 4. 5. 6. 7. 8. 9. 10.

500

B. S. Rinkevicius, Laser Diagnostics of Flows [in Russian], Izd. MEI, Moscow (1990). V. K. Andreev et al., Modern Mathematical Models of Convection [in Russian], Fizmatgiz, Moscow (2008). O. A. Evtikhieva, I. L. Raskovskaya, and B.S. Rinkevicius, Laser Refractography [in Russian], Fizmatgiz, Moscow (2008). Yu. A. Kravtsov and Yu. I. Orlov, Geometric Optics of Inhomogeneous Media [in Russian], Nauka, Moscow (1980). O. A. Evtikhieva, Izmer. Tekhn., No. 5, 35 (2006); Measur. Techn., 49, No. 5, 446 (2006). K. M. Lapitskii, I. L. Raskovskaya, and B. S. Rinkevicius, Izmer. Tekhn., No. 7, 28 (2008); Measur. Techn., 51, No. 7, 740 (2008). K. L. Lapitskii, Metrologiya, No. 7, 33 (2008). M. B. Vinogradova, O. V. Rudenko, and A. P. Sukhorukov, Theory of Waves [in Russian], Nauka, Moscow (1990). V. I. Artemov et al., Proc. 4th Russian National Conf. On Heat Exchange, Izd. Dom MEI, Moscow (2006), Vol. 3, p. 42. K. M. Lapitskii, I. L. Raskovskaya, and B. S. Rinkevicius, Abstr. of 12th Intern. Symp. on Flow Visualization (J. Kompenhans, A. Schroeder, and I. Grant, eds.), German Aerospace Center (DLR), Goettingen (2006), p. 24.