Äîêóìåíò âçÿò èç êýøà ïîèñêîâîé ìàøèíû. Àäðåñ îðèãèíàëüíîãî äîêóìåíòà : http://lizard.phys.msu.su/home/science/Dmitriev-Ivanov-Xoxlov-11-JMathSci-Diffuser.pdf Äàòà èçìåíåíèÿ: Thu Feb 10 10:31:27 2011 Äàòà èíäåêñèðîâàíèÿ: Mon Oct 1 20:07:33 2012 Êîäèðîâêà:

Journal of Mathematical Sciences, Vol. 172, No. 6, 2011

NUMERICAL SIMULATION OF LIGHT PROPAGATION THROUGH A DIFFUSER A. V. Dmitriev, A. V. Ivanov, and A. R. Khokhlov UDC 538.935+517.958+537.84

Abstract. Diffusers are important elements of many illumination systems, for example, in computer and mobile phone displays or advertising panels, etc. In this article, the light propagation in a diffuser with optically soft inclusions is described with the help of the Fokker­Planck equation, i.e., a transfer equation with a diffusion term in the space of radiation-propagation directions. The coefficient of angle diffusion is calculated using the Mie theory. The equation is solved numerically using the stochastic analog method, and the space and angle distribution of the radiation that passed through the diffuser is calculated. The results can be useful for optimization of diffuser parameters, and the method can be applied to many problems of turbid media with optically soft particles.

Very often, point or linear light sources are used for the illumination of large transparent pictures and special devices are used to ensure a uniform illumination. Photo, cinema, and other pro jectors usually contain special lens condensers for this purpose, but in such thin and compact devices as LCDs of TV sets, computers, and mobile phones, so-called diffusers are used instead of lens condensers. A diffuser typically consists of one or a few layers of transparent glass or plastic, each containing microscopic particles made of other material with different sizes and in different concentrations, which scatter the light. The particle diameters are usually of order of the light wavelength. Similar diffusers can be found also in advertising windows, where they are made routinely of a single sheet of white plexiglas placed between luminescent lamps and an image to smoothen illumination nonuniformity. Although it is a rather simple device, a diffuser has to fulfill many conditions that may contradict to each other or be completely different for different devices. To ensure a highly uniform illumination, it should be optically dense, but to avoid brightness reduction, it should be transparent. In advertising and TV sets, the diffusers must ensure as wide an angle of observation as possible, but the displays of ATMs and business notebooks must be visible only to the person situated in front of them. Thus, methods of creating diffusers with different properties are necessary and hence the problem of modeling light distribution in diffusers is important. However, from a scientific point of view, this problem belongs to a very complex part of optics: the optics of turbid media. A turbid or scattering medium is a uniform and transparent substance (the matrix) that contains particles of another substance, and these particles scatter light propagating in the media. The simplest turbid media is fog, i.e., small water droplets in air. A diffuser is also a turbid medium. The theory of light propagation in turbid media is based on the so-called radiative transfer equation [1, 4, 10, 12], but methods for its analytical solution are lacking, and it is also not easily solved numerically. For our present computational purposes, taking into account a nonuniform light-intensity distribution in three-dimensional space, it is more convenient to write the radiative transfer equation in the form [12] 1 1 I (t, n, r) + nI (t, n, r) = cm t lsc 1 d x( )I (t, n , r) - I (t, n, r). 4 lex (1)

Here I is the light intensity, t is the time, r is the coordinate vector, cm is the speed of light in the diffuser media, n and n are the unit vectors in the light-propagation directions, integration is performed over all n directions, x( ) is the light-scattering indicatrix (the probability of light scattering in a given direction),
Translated from Fundamentalnaya i Prikladnaya Matematika, Vol. 15, No. 6, pp. 33­41, 2009.

782

1072­3374/11/1726­0782 c 2011 Springer Science+Business Media, Inc.


being the scattering angle (the angle between n and n ), lsc and l respectively: 1 1 , lex = , lsc = Nsc Nex

ex

are scattering and extinction lengths,

N is the concentration of scatterers (inclusions in the diffuser media), and sc and ex are scattering and extinction cross sections of a single scatterer. The light frequency enters the equation as a parameter that influences the shape of the scattering indicatrix as well as the values of the light velocity in the diffuser media and the cross sections. The radiative transfer equation is a complex integro-differential equation in a 6-dimensional hyper-space, its components being time, two components of n, and three components of r, or in a 5-dimensional hyper-space if the stationary case is considered, and effective methods of its solution in the spatially nonuniform case are lacking. But if the refractive indices of the media and of the scatterers are close to each other (so-called "optically soft scatterers"), which is precisely the case for most production diffusers, then the equation can be simplified with the help of the method developed by V. V. Sobolev [11, 12]. This method is often called the small-angle diffusion approximation. The method is based on the fact that for an optically soft scatterer the light scattering is relatively ineffective and hence the scattering indicatrix is a narrow peak around the incident-light-ray direction. The more optically soft the scatterer, the narrower the indicatrix, and in the limiting case where the difference of the refractive indices becomes infinitely small, the indicatrix tends to the angle delta-function. This means that for optically soft scatterers the mean change of the light-propagation direction in an elementary scattering event is small and hence complete angle relaxation of the incident ray requires many scattering events. This is a feature of the so-called slow processes [5], and they can be described by the differential Fokker­Planck equation instead of more complicated integral equations [9]. For the radiative transfer equation, the corresponding transformation can be performed in the following way [11, 12]. If the scattering indicatrix is narrow, the difference between n and n directions is small on the average, and hence the difference between I (n) and I (n ) is small, so that one can substitute instead of I (n ) its Taylor expansion in terms of - and - , where and are the direction angles of the vector n, and and are those of n . As is usual in the Fokker­Planck theory, it is sufficient to keep the expansion terms up to the second order and neglect terms of third and higher orders [5, 9]. Substituting the resulting expansion into the integral on the right-hand side of the radiative transfer equation, one obtains 1 I 1 2I d . x( )I (t, n , r) = I (t, n, r)+ D sin + 4 sin sin2 2 The first term on the right-hand side comes from the zeroth term of the Taylor expansion, where the condition d x=1 4 has been taken into account, the first-order terms disappear after the integration, and the term with the second derivatives of I comes from the second-order terms of the expansion, the dimensionless coefficient D being the second angle moment of the scattering indicatrix: 1 D= 8


d x( )sin3 .
0

The term with derivatives describes diffusive motion of the light-propagation direction in the angle space, and D has the sense of the angle diffusion coefficient of light in the diffuser due to the scattering by inclusions. It was shown in [11, 12] that the angle diffusion coefficient can be expressed also via x1 , the 783


first coefficient of the scattering indicatrix expansion in terms of the Legendre polynomials: 1 D = (3 - x1 ), 6 3 x1 = 2


d x( )cos sin .
0

We use the latter equalities for calculations of D. Using the well-known expressions of differential operators in the spherical coordinate system [8], one can rewrite the second-order differential term on the right-hand side of the radiative transfer equation in a more compact vector form: I 1 sin sin + 1 2I = - div([n â [n â grad]]I ) = -([n â [n â]]I ), sin2 2

which allows us, finally, to write down the radiative transfer equation in the shape of the Fokker­Planck equation in the angle space: 1 I (t, n, r) + nr I (t, n, r)+ cm t 1 l
ex

-

1 lsc

I (t, n, r) = -

D n {[n â [n ân ]]I (t, n, r)}. lsc

This is now a differential equation instead of the original integro-differential one, and it can be effectively solved numerically using the stochastic analog method, a variant of the Monte Carlo method. The equation is valid inside the layer bounded by planes rz = 0 and rz = Lz . At the layer boundaries, Snell's law and total internal reflection and partial reflection, described by Fresnel formulas, are taken into account. Outside the layer, geometrical optics is used. The inclusions in the diffuser medium are assumed to be spherical in shape, and their size distribution is assumed to be a gamma-distribution [2, 12]. The corresponding scattering indicatrix, which is used for calculation of the angle diffusion coefficient, is calculated using the Mie theory [3, 6, 13]. We are interested in the time-independent situation. The steady-state solution of the radiative transfer equation is found by using the Monte Carlo method (the stochastic analog method) [7, 14]. Briefly, rays with random directions that originate from the light source and then spread through the system, experiencing scattering and reflection with proper probabilities, are chosen. Their exit points from the diffuser are registered together with their exit directions and used for tracing the lateral and angle intensity distribution, the intensity being proportional to the lateral distribution function of rays at the exit surface of the diffuser. The exit angles of the rays are used to obtain the angle structure of the intensity of the radiation that passed through the diffuser with spatial resolution. It is worth noting that light reflection and refraction as well as total internal reflection are easily accounted for in this method, which distinguishes it favorably from the well-known doubling method. Another advantage is that the method naturally gives space-resolved angle distributions. The numeric scheme for one ray is as follows: rk nk
+1 +1

= rk + hcnk , = R R(nk , 23 k ) ik â nk , 2 -2hD ln 3k |ik â nk |
+1

cos 23k

+2

nk ,

where k is the number of the computation layer, h is the computation step, ik is the unit vector, corresponding to the direction of the minimum component of the vector nk , {j } is a sequence of random digits uniformly distributed on the interval (0, 1), and R(o, ) is the rotation matrix of the unit vector o by angle according to the corkscrew rule: 2 2sin 2 ox oz sin 2 + oy cos 2 (ox - o2 - o2 )sin2 2 +cos2 2 2sin 2 ox oy sin 2 - oz cos 2 y z R(o, ) = 2sin 2 ox oy sin 2 + oz cos 2 (o2 - o2 - o2 )sin2 2 +cos2 2 2sin 2 oy oz sin 2 - ox cos 2 . y x z 2sin 2 oy oz sin 2 + ox cos 2 (o2 - o2 - o2 )sin2 2 +cos2 2 2sin 2 ox oz sin 2 - oy cos 2 z x y 784


To calculate the distribution function at the exit surface of the diffuser having area [-Lx /2,Lx /2] â [-Ly /2,Ly /2], a rectangular uniform mesh with Nx â Ny cells is used. An isotropic mesh that is an inscribed polytope of a unit sphere is used in every cell to trace the angle distribution function. Each face is close to a regular-shaped triangle. The basis for mesh generation is the dodecahedron. Each face of the dodecahedron is divided into five equal close-to-regular-shaped triangles (60 triangles in total). If needed, the mesh can be further divided to the requested degree by recursive fragmentation of each triangle into four parts. There are 60 · 4r cells in the sphere mesh, r being the number of divisions. The angle distribution function can be displayed with the help of a special viewer y = 1 - 2/ , x = (1 - / )(1 - y 2 ), transformation, where and are the spherical coordinate angles and x and y are coordinates in the plane. This pro jection is a certain prototype of those used in designing world maps. Furthermore, the averaged-over-the-azimuthal-angle direction diagram is calculated as a function 2 of on the x2 + y 2 < r section of the diffuser and on the uniform mesh [0, ] with N cells. Also, the diffuser transmittance factor Ieff , which indicates the portion of radiation that has passed through the diffuser, is calculated. In order to get some acceptable results it is necessary to obtain sufficiently large statistics. The number of rays N used in the calculations should be much greater than the number of cells used for tracing the distribution function: N Nx Ny · 60 · 4r . It is worth noting that the longer the ray, the less its intensity because of permanent scattering and reflections, so it is useless to follow rays that are too long; in our calculations, we limited the ray length to lmax (20­100)Lz . The numeric scheme is realized as a high-performance software system, written in C++, Python, and Fortran languages for the Linux OS. The system is supplied with a user-friendly interface. This interface provides startup, continuation or "cloning" of calculations, graphic visualization of results, and automatic export of figures into generally accepted .eps and .gif graphics formats. Calculation data are also stored in files in a special format. Figures 1­4 show the results of the light-propagation modeling through a plexiglas diffuser with microscopic particles of silica glass.

Fig. 1. The radiation intensity that comes to an observer from the surface of the diffuser illuminated by a linear light source under its center (the bright bar below).

785


Fig. 2. The radiation intensity that comes to an observer from the surface of the diffuser illuminated by a point light source from below (the source is not visible).

Fig. 3. The angular radiation distribution in the center of the diffuser illuminated by a point light source under the diffuser center. At the top, the intensity of the light that passed through the diffuser is seen; at the bottom, the radiation reflected from the diffuser lower boundary.

Fig. 4. The angular radiation distribution at the edge of the diffuser illuminated by a point light source under the diffuser center. At the top: the light that passed through the diffuser; at the bottom: the reflected light.

786


The method developed here and the corresponding program complex can be useful for calculating optimal diffuser parameters for a particular device as well as for modeling light transfer in different kinds of turbid media with boundaries of complex geometry, rough boundaries, etc. Acknowledgment. This work was supported by the Federal special-purpose program "Scientific and scientific-educational personnel of innovative Russia in the years 2009­2013" via contract P-2312. REFERENCES 1. L. A. Apresyan and Yu. A. Kravtsov, Theory of Radiative Transfer. Statistical and Wave Aspects [in Russian], Nauka, Moscow (1983). 2. V. A. Babenko, L. G. Astafyeva, and V. N. Kuzmin, Electromagnetic Scattering in Disperse Media, Springer-Verlag, Berlin (2003). 3. C. F. Bohren and D. R. Huffman, Absorption and Scattering of Light by Smal l Particles, Wiley, New York (1983). 4. S. Chandrasekhar, Radiative Transfer, Dover, New York (1950). 5. A. V. Dmitriev, Basics of the Statistical Physics of Materials [in Russian], Nauka, Izd. Mosk. Univ., Moscow (2004). 6. H. C. Van de Hulst, Light Scattering by Smal l Particles, Wiley, New York (1957). 7. A. V. Ivanov, "Kinetic modelling of the dynamics of magnetics," Mat. Model., 19, No. 11, 89­104 (2007). 8. G. A. Korn and T. M. Korn, Mathematical Handbook, McGraw-Hill, New York (1968). 9. E. M. Lifshits and L. P. Pitaevskii, Physical Kinetics, Pergamon, Oxford (1981). 10. M. M. Mishchenko, L. D. Davis, and A. A. Lacis, Scattering, Absorption, and Emission of Light by Smal l Particles, Cambridge University Press (2002). 11. V. V. Sobolev, Dokl. Akad. Nauk SSSR, 177, No. 4, 812­815 (1967). 12. V. V. Sobolev, Light Scattering in Planet Atmospheres [in Russian], Nauka, Moscow (1972). 13. W. J. Wiscombe, "Improved Mie scattering algorithms," Appl. Optim., 19, 1505­1509 (1980). 14. G. I. Zmievskaya, "Numerical stochastic models of nonequilibrium processes," Mat. Model., 8, No. 11, 3­40 (1996). A. V. Dmitriev, A. V. Ivanov, and A. R. Khokhlov Moscow State University, Moscow, Russia E-mail: dmitriev@lt.phys.msu.su

787