## Abstract

Modern optical diagnostics for quantitative characterization of polydisperse sprays and other aerosols which contain a wide range of droplet size encounter difficulties in the dense regions due to the multiple scattering of laser radiation with the surrounding droplets. The accuracy and efficiency of optical measurements can only be improved if the radiative transfer within such polydisperse turbid media is understood. A novel Monte Carlo code has been developed for modeling of optical radiation propagation in inhomogeneous polydisperse scattering media with typical drop size ranging from 2 *μm* to 200 *μm* in diameter. We show how strong variations of both particle size distribution and particle concentration within a 3D scattering medium can be taken into account via the Monte Carlo approach. A new approximation which reduces ~20 times the computational memory space required to determine the phase function is described. The approximation is verified by considering four log-normal drop size distributions. It is found valid for particle sizes in the range of 10–200 *μm* with increasing errors, due to additional photons scattered at large angles, as the number of particles below than 10 *μm* increases. The technique is applied to the simulation of typical planar Mie imaging of a hollow cone spray. Simulated and experimental images are compared and shown to agree well. The code has application in developing and testing new optical diagnostics for complex scattering media such as dense sprays.

©2005 Optical Society of America

## 1. Introduction

Non-intrusive optical diagnostics of inhomogeneous turbid media (e.g. sprays, aerosols, smoke, fog, foams, *etc*) are of interest in many research domains including combustion engineering, meteorology, biomedicine, *etc*. The major issues occurring during these optical measurements are frequently related to the multiple scattering of the initial laser radiation with the surrounding scattering particles. For instance, the accuracy in optical measurements of a fuel spray depends on both the droplets concentration and on the dimension of the given sprays. These two parameters are directly responsible on the amount of multiple scattering that will occur [1]. At low droplet concentration and small dimension of the scattering medium, it can be assumed that the probe beam retains essentially the same intensity as it traverses the spray, and that each photons packet scatters from only one droplet. Under these conditions interpretation of detected signal is straightforward. However at high droplet concentrations and large dimensions of the scattering medium a number of effects cause errors if uncorrected:

- The probe beam is attenuated as it traverses the spray [2]. Depending on their position in the spray, droplets are not all lighted with the same initial intensity.
- The scattered light is also attenuated by “secondary scattering” from droplet lying between the probe beam and the detector.
- Some “extraneous light” is also detected after being multiply scattered by the surrounding droplets.

Finally, a part of the detected radiation have scattered from more than one droplet and carry information about all the droplets they have encountered. The magnitude of each error varies with position in the image, in a manner dependent on the spray structure. Correction is not simple as of course the structure of the spray is unknown.

However, the accuracy and the efficiency of modern optical diagnostics can be predicted and improved if the radiative transfer within these turbid media is understood [3]. The Monte Carlo (MC) method is one of the most versatile and well developed probabilistic techniques used for the simulation of optical radiation propagation through various complex scattering media, such as the atmosphere [4, 5], biological tissues [6], geological structures [7], and sprays [8]. Nowadays MC modeling is widely applied to optimize the source-detector configuration of a particular experimental set-up [9 ,10], to estimate fluence rate distribution within the human tissues [11], to predict detected optical signal/reflectance spectra [12], to analyze coherent effects in multiple scattering [13] and for other theoretical studies. Most of existing MC models [3–13] considers monodisperse homogeneous scattering media, or media combining several monodisperse homogeneous layers. In practical sprays density and size distribution of droplets vary strongly within the medium and the assumption of a homogeneous medium cannot be applied. In this paper a unique Monte Carlo code which deals with polydisperse and highly inhomogeneous media is presented. We show and verify a technique for determining the local scattering phase function, able to take into account the Mie phase function of a wide range of diameters, with low computational requirements. This technique finds its application in the study of highly inhomogeneous polydisperse turbid media. In this most complex scattering case, the use of only one phase function (representative of a single particle only or averaged over a particle size distribution) cannot be performed.

The code is applied to the simulation of a typical laser sheet imaging of industrial hollow cone spray. Finally, the use of the code for improving and developing new optical diagnostics of dense sprays is discussed.

## 2. Modeling of light propagation in inhomogeneous and polydisperse scattering media

The main assumption of the MC technique is to define the light emitted by the source as point entities (photon packets), here called photons for ease. Each photon enters the medium which contains scattering and absorbing centers (droplets or particles) with an initial direction and is tracked as it travels through the medium. The trajectory of the photons is governed by probability density functions beforehand defined: The probability that a photon is scattered, the probability that it is absorbed and the probability to follow a new direction of propagation after a scattering event. The principle steps of the MC technique are the followings: The free path length *l* before each light-particle interaction is derived from the Beer-Lambert law and is calculated as a function of the extinction coefficient *μ _{ext}* using a random number ξ uniformly distributed between 0 and 1:

*l*= - ln(

*ξ*)/

*μ*with

_{ext}*μ*=

_{ext}*N*

^{*}

*σ*.

_{ext}*N*is the number density of scattering particles and

*σ*is the extinction cross section. At each interaction with a particle photons can be either absorbed or scattered depending on the medium albedo Λ. If particles are non-absorbing, the extinction coefficient is then equal to the scattering coefficient(

_{ext}*μ*=

_{ext}*μ*+

_{sca}*μ*where

_{abs}*μ*is the scattering coefficient and

_{sca}*μ*is the absorption coefficient) and the albedo equal 1 (Λ =

_{abs}*μ*/(

_{sca}*μ*+

_{sca}*μ*)). Other numerical methods for the determination of absorption and scattering of optical radiation by particles have been described and discussed in [16]. In the MC technique, the scattering phenomena are assumed independent one to each other. This requires a distance between particles of greater than three times the radius [17]. After a scattering event, the photons new direction is selected with a random number and a Cumulative Probability Density Function (CPDF) calculated from the appropriate scattering phase function

_{abs}*f*. The scattering phase function defines the relative light intensity scattered in all directions as a function of the incident wavelength, the refractive index of the surrounding medium, the state of polarization and the properties of the particle encountered (size, shape, orientation, constitution and refractive index). Depending on the particle size parameter, Mie [17, 18], Rayleigh-Gans [17, 18] or Henyey-Greenstein [19] phase functions are typically used. The polar scattering angle

*θ*defined between 0 and

_{s}*π*is found from the inverse CPDF of

*f*by

*θ*= CPDF

_{s}^{-1}(ξ) (where ξ is a random number between 0 and 1). The azimuthal scattering angle

*φ*is uniformly distributed between 0 and 2

_{s}*π*. When a new direction of propagation is defined, the position of the next scattering point is calculated again and the process is repeated until the photon is either absorbed or leaves the medium. The total number of photons sent depends on the accuracy desired and on the characteristics of the detection. The final direction of propagation, the final position, the number of scatters, and the total path length are calculated at the end of each photons journey. If the conditions of detection are met (e.g. photon lies within the field of view of the detector with its trajectory within the acceptance angle), these data are written to disk. The process is repeated for a large amount of photons such that the distribution of all light intensity impinging on the detector has been found in the 3D coordinate system. If an infinite number of photons were sent, the exact solution of the Radiative Transfer Equation (RTE) would be obtained. The RTE is the mathematical expression which describes the conservation of radiant energy of optical radiation through a turbid media [16]. For most practical geometries the RTE cannot be solved analytically. The MC technique can handle all conceivable geometrical configurations of source, medium and detector and is the most flexible method of reaching an accurate approximate solution. Thus the MC technique is a powerful tool to understand and simulate scattering processes in different turbid media.

The complexity of the structure of a turbid medium has a direct bearing on the complexity of the MC model required for the simulation. Turbid media are characterized by the concentration of particles, their distribution in space, and by the number of different particle sizes or types present. Five cases of turbid media are identified Table 1 depending on the homogeneity of particle size/type and on their spatial distribution. The simplest case corresponds to homogeneous and monodisperse media: As only one type of particle is considered, only a single scattering CPDF and a single extinction cross section *σ _{ext}* are required. If a homogeneous medium is polydisperse, the scattering CPDF used must be deduced from the phase function averaged over the total distribution of particle sizes/types. The extinction coefficient

*μ*is in this second case calculated from the average extinction cross section

_{ext}*σ*¯

_{ext}. By definition, homogeneous media are defined by a constant number density

*N*of particles in every single point of the sample. When the number density of particles varies from place to place the medium becomes inhomogeneous. Both the extinction coefficient and the scattering CPDF change with location. Working with an inhomogeneous medium requires the scattering medium to be decomposed into elementary volumes in which these properties are homogenous. In the presented model these elementary volumes are cubic cells of constant size (Fig.1). The size and the number of the cells are chosen based on the accuracy required and on the geometry of the medium. The path length between scattering events of a photon transferring from one cell to another is corrected proportionally to the ratio between the extinction coefficients of the “last cell” crossed and the extinction coefficient of the “new cell” encountered. If

*μ*

_{ext}(new cell) <

*μ*

_{ext}(last cell) , the free path-length

*l*is increased, and if on the contrary

*μ*

_{ext}(new cell) >

*μ*

_{ext}(last cell) ,

*l*is decreased.

Inhomogeneous media can be monodisperse, uniformly polydisperse or polydisperse. For monodisperse (only one size of particle present) and uniformly polydisperse (particle size distribution constant with location) inhomogeneous media, the scattering process and hence the related CPDF is assumed identical in every cell. The variation of the extinction coefficient *μ _{ext}* is simply related to the variation of the number density of particles.

For inhomogeneous polydisperse media, both number density of particles and particle size distribution varies with location. The extinction cross section and the scattering CPDF must therefore to be defined in each cell (Fig. 1). This constitutes the most complex scattering case. Introducing a scattering CPDF for each cell requires a large amount of input data, particularly if the medium is represented with a large number of cells. Generating the correct scattering CPDF at each scattering event dramatically increases the running time. One solution to this problem uses several averaged scattering CPDF (stored in lookup tables), each one representing either the scattering of a set of particles of similar size, or representing a particular size distribution. This approximation is verified against the more rigorous approach in the next section.

## 3. Verification of the phase function approximation

#### 3.1. Calculation methods

Four log-normal distributions of particle sizes have been chosen using different values of average diameter *d*¯ and standard deviation *σ*. The log-normal function is defined by:

Where *T* is the shape parameter and *M* is the scale parameter. The mean diameter and the standard deviation are respectively given by:

The distributions tested are based on mean diameters *d*¯ = 5 *μm* and *d*¯ = 40 *μm* with standard deviations *σ* equal to 10% and 80% of *d*¯ (Fig. 2). The bin width is 0.4 /*μm*. Distribution (a) is representative of an automotive fuel injector spray, distribution (c) of a medical nebulizer spray, and the other distributions are included to show the effect of standard deviation.

Two methods for representing the local scattering phase function in polydisperse homogeneous media were tested. The method one (M_{1}) is based on the determination of the
average phase function *f*¯ over the complete distribution of drops size Eq. (4).

Where *n*(*d*) is the number of drops of diameter *d*. *f*¯ has been calculated for the four
particle distributions described above and the scattering CPDF of *f*¯ is deduced for each of these distributions. Only one scattering CPDF representative of the complete drops distribution is used in the MC simulations with the method M_{1}. In the case of infinitesimal size bin width of the dropsize distribution M_{1} would give the exact mathematical solution of the global scattering process that occurs in a homogeneous polydisperse turbid media with independent scattering.

In method two (M_{2}), 25 different scattering CPDF are defined such as each CPDF is representative of the scattering by particles belonging to a class of drop sizes. Even if several approaches can be employed to determine these CPDF, only one has been found valid. The first approach consists to takes into account only the phase function corresponding to the middle drop size of the class bin. As the scattering phase functions do not change linearly with the drops diameter this method has been rejected. The second approach consist to add the phase function of the minimum value of the bin to the one corresponding to the maximum value and to divide the result by 2. Once again this approximation has been rejected due to the non-linear changes of the phases function with droplet sizes. Finally the appropriate method is based on calculation of the phase function averaged over the range of particle size, with an equal number of drops for each given size but weighted with the corresponding scattering cross section (Eq. (4) with *n*(*d*) equals 1 for every *d*). When a scattering event occurs, the diameter of the particle encountered is determined with the distribution of droplet size and a random number. The probability *P*(*d*
_{1}) of a drop of diameter *d*
_{1} being encountered by a photon packet is given by:

Once the diameter of the particle encountered is found, the correct approximate scattering CPDF can be chosen. Note that if the exact phase function (corresponding to the drops size reached) was chosen (instead of the approximated one), method M_{2} would be equal to method M_{1}. The accuracy of M_{2} is then directly related to the difference between the real scattering CPDF of the droplet encountered and the approximate CPDF chosen. Reducing the size range over which scattering CPDF are averaged increases the accuracy of the technique. In the present work 25 classes of drop size are used and the range of sizes for each class is varied according to the rate of change of scattering CPDF with droplet size. The boundaries between each size class are selected by hand to minimize the difference between the CPDF on neighboring classes.

It is seen in Fig. 3 that for small particles neighboring scattering CPDF diverge strongly and overlap at small angles *θ _{s}* when

*d*≤ 15

*μ*m. For large particles, the scattering CPDF do not overlap and remain close one to each other even when the particle size interval is large. Thus to maximize the accuracy of the technique while minimizing the memory requirements, the range of sizes is kept small for small droplets (~1

*μ*m) and large (up to 29

*μ*m) for large droplets. Further reducing the particle size increases the accuracy of the technique at the expense of greater memory requirements.

#### 3.2. Description of the simulation

Each droplet size distribution in Fig. 2 has been used. The droplets are contained in a homogeneous single cubic cell of dimension *L* = 50 *mm*. A cylindrical flat laser beam *S* of 20 *mm* diameter enters through the scattering sample crossing perpendicularly the Y=0 plane (back face) and exiting through the Y= *L* plane (front face) (Fig. 4).

The source wavelength *λ* is 532 *nm* and the light is assumed unpolarized. Intensity profiles on the back face (back scattering) and on the front face (forward scattering) are recorded for different detector acceptance angles *θ _{a}*. In each simulation 100 millions photons are sent. The surrounding medium is air (refractive index equals 1+0.0

*i*). The droplets are spherical and non-absorbing with refractive index 1.4+0.0

*i*and the scattering CPDF are calculated from the Mie theory. Scattering and extinction coefficients are then equal and the simulations are run with

*μ*fixed to 0.12 and 0.24

_{sca}*mm*

^{-1}. The resulting optical depths are respectively 6 and 12, corresponding to the intermediate single-to-multiple scattering regime. The range of particle size is 2

*μ*m − 200

*μ*m (typical of droplet sizes in fuel sprays) and the resulting size parameter

*S*

_{p}=

*λ*

^{*}

*π*/

*d*is: 11.81 ≤

*S*

_{p}𢙔 1181.05.

#### 3.3. Results and comparison

Results obtained from M_{2} are compared to the results obtained from M_{1} by analyzing quantitatively the light intensity distribution on the front and back face of the scattering medium (Fig. 4). Note the intensity scale is different for each image. The images presented are obtained using M_{1} for a single dropsize distribution with an average diameter of 40 *μm* and a standard deviation of 32 *μm* (Fig. 2(a)). The scattering coefficient *μ _{sca}* is fixed to=0.12

*mm*

^{-1}and 0.24

*mm*

^{-1}giving respectively an average number of scatters per photon of ~6 and ~12 scatters. Figure 5(a) (

*μ*=0.12

_{sca}*mm*

^{-1}) and (b) (

*μ*=0.24

_{sca}*mm*

^{-1}) demonstrate the broadening of the beam on the front face as the optical depth increases. The images illustrate the quantitative distribution of the forward and backward scattered light considering all scattering orders (Fig. 5(a), (b), (e) and (f)) and with the detector filtered to detect single scattering only (Fig. 5(c), (d), (g) and (h)). Figure 5(a) and (b) show that doubling

*μ*the forward light intensity is strongly attenuated (by a factor of ~3.3) and that the shape of the laser beam is no longer clearly defined. On the contrary for back scattering (Fig. 5(e) and (f)), as

_{sca}*μ*increases the detected intensity increases also but the pattern of scattered radiation does not change significantly.

_{sca}Single scattering detected in the forward direction shows a faithful reconstruction of the laser beam for *μ _{sca}* = 0.12

*mm*

^{-1}(Fig. 5(c)). However the intensity of single scattering is weak compared to the amount of multiple scattering for

*μ*= 0.24

_{sca}*mm*

^{-1}(Fig. 5(d)). It can be seen from Fig. 5(g) and (h) that single back scattered signal remains relatively constant for both scattering coefficients. The effects of the detection acceptance angle are also investigated. In Fig. 5, all photons reaching the detection areas are detected (acceptance angle

*θ*= 90°). As found in other simulations [9], the acceptance angle can be used to optimize the ratio of singly to multiply scattered photons detected.

_{a}In Fig. 6 the detection acceptance angle *θ _{a}* is reduced to 5°.

*μ*=0.12

_{sca}*mm*

^{-1}in Fig. 6(a) and (c), corresponding to Fig. 5(a) and 5(c), and 0.24

*mm*

^{-1}in Fig. 6(b) and 6(d), corresponding to Fig. 5(b) and 5(d). The total intensity on the front face is strongly reduced but the boundaries of the laser beam appear clearly (Fig. 6(a) and (b)). The single scattering intensity detected per pixel remains close for both acceptance angles

*θ*= 90° Fig. 5(c) and 5(d) and

_{a}*θ*= 5° Fig. 6(c) and (d). This indicates that most of single scattered photons propagate with a polar scattering angle

_{a}*θ*less than 5° and shows the high intensity of scattering in the forward direction in Mie scattering processes.

_{s}The scattering coefficient, the geometry of the sample and the scattering phase function all influence the number of scattering events *n* occurring and the total path length *L* of the photon packets. Slice differences in the scattering phase functions (between the exact one and the approximated one used in the simulations M_{2}) do not affect significantly the parameters *n* and *L*, but can modify the final intensity distribution. The images presented in Fig. 5 and Fig. 6 have been obtained from method M_{1}. Corresponding images have been generated by applying M_{2}. Due to the symmetry of the images, only the intensity profile along a line passing from the centre of the laser beam X=25 *mm* until the edge of the image X=50 *mm*, is considered with Z fixed to 25 *mm* at Y=0 or Y=*L*. The comparison between the two methods is made by calculating the ratio of intensities along this line of the image generated with M_{1} to the image generated with M_{2}. This ratio is plotted (Fig.7) for *μ _{sca}*= 0.12

*mm*

^{-1}and

*μ*= 0.24

_{sca}*mm*

^{-1}using once again the log-normal droplet size distribution defined by

*d*¯ = 40

*μm*and σ = 32

*μm*(Fig. 2(a)). Total intensity (single and multiple scattering taken together) are detected on the front face and the back face of the scattering cube with a detection acceptance angle

*θ*= 90°. It is seen from Fig.7 that the ratio M

_{a}_{1}/M

_{2}(Image method one / Image method two) remains equal to ~1 (±0.02) for both forward and back light scattering. These results show the very good agreement between the two methods at large detection acceptance angles.

Figure 8 shows the same data, but this time with a detection acceptance angle *θ _{a}* of 5°. In Fig. 6 it is seen that with this restricted acceptance angle the number of detected photons is very low outside of the projected area of the incoming beam. This low photon count makes the M

_{1}/M

_{2}intensity ratio noisy in Fig. 8 at X > 35

*mm*. At X < 35

*mm*the greater number of detected photons give a better defined ratio. Here it is seen that the results from M

_{2}match the results from M

_{1}for small detection acceptance angles as well as for large. Doubling the scattering coefficient from 0.12 mm

^{-1}to 0.24 mm

^{-1}significantly reduces the number of photons reaching the front face and making the M

_{1}/M

_{2}intensity ratio noisier for

*μ*= 0.24

_{sca}*mm*

^{-1}Fig. 8(b) than for

*μ*= 0.12

_{sca}*mm*

^{-1}Fig. 8(a).

In Fig. 9 the comparison is performed for single scattering detection with *μ _{sca}*= 0.12

*mm*

^{-1}. It is seen that for both acceptance angles 5° and 90° the ratio fluctuates equally either side of 1 when X < 35

*mm*. At large distance from the laser beam centre (X > 35

*mm*) results diverge between the two detection apertures: If

*θ*=90° (Fig. 9(a)), few photons are detected giving strong statistical fluctuations in the resulting ratio M

_{a}_{1}/M

_{2}. If however

*θ*=5° (Fig. 9(b)), no singly scattered photons are detected for either method and a flat line is plotted for X > 40

_{a}*mm*. A small number of photons are detected at the edge of the laser beam at the intermediate distance 35

*mm*< X < 40

*mm*giving a noisy interval.

These comparisons show that the differences between the two methods are only observed when the detected signal is weak. These differences are caused by the strong statistical fluctuations which occur when the amount of collected data resulting from probability laws is too low. However where these fluctuations are strong, the ratio is biased to values greater than one. More photons are then detected with M_{1} than M_{2} on the front face when the signal is weak. It is deduced that the weight given to the scattering phase function of large particles (with large forward scattering lob) is then more important in M_{1} than in M_{2}.

These results demonstrate that apart from very small differences in the number of detected photons scattered at high angles when the signal is weak, the results obtained with the phase function approximation used in method M_{2} are in excellent agreement with the rigorous method M_{1}, for droplet sizes over the range 2–200 *μm* of particle size, for different *θ _{a}*=5° and 90° and with

*μ*varying over a factor of 2.

_{sca}It remains to verify method two for other particle size distributions. Fig. 10 shows the results for the three other log-normal distributions, defined respectively by *d*¯ = 40 *μm* with σ = 4 *μm* (Fig. 10(a)), *d*¯ = 5 *μm* with σ = 4 *μm* (Fig. 10(b)) and *d*¯ = 5 *μm* with σ = 0.5 *μm* (Fig. 10(c)). Photons from all scattering orders have been detected with an acceptance angle *θ _{a}* = 90° on the front face, assuming a scattering coefficient of 0.12

*mm*

^{-1}.

When the average diameter is 40 *μm* with 4 *μm* standard deviation (σ = 10% of *d*¯) Fig. 10(a)) the ratio of intensities from the two methods is once again ~1 (±0.02). This result verifies the use of method two for distributions of particles based on large mean diameters with small relative standard deviation. Greater differences between method two and method one appear when small drops are considered. It is seen that the ratio M_{1}/M_{2} is greater than 1 and reaches a maximum of ~1.2 when *d*¯ = 5*μm* with σ = 4 *μm* (Fig. 10(b)). If the standard deviation is reduced to 0.5 *μm* (Fig. 10(c)) the differences between the two methods are increased with a ratio lying between 1.15 < M_{1}/M_{2} < 1.4. It is deduced from this results that method one gives more forward light scattering than method two when only small particles are considered. These results were expected from Fig. 3 due to the differences in the averaged scattering CPDF for small drop sizes which result from the averaging used in the two methods. Referring to the size distributions plotted in Fig. 2, the comparison above demonstrates that the phase function approximation used in M_{2} gives accurate results for any distribution of spherical drops comprised between 10 and 200 microns, except for very small differences in the number of photons scattered at high angles where the singly scattered signal is weak. For particles smaller than 10 *μm*, discrepancies in the global light intensity distribution appear between the exact and the approximated solution. These differences can be corrected by reducing the particle size interval used for small drops (*d* < 10 *μm*) in M_{2} (Fig. 3). Particles from 2 to 200 *μm* have been considered with a step of 0.4 microns. The exact approach would require using a total of ~500 phase functions. The use of only 25 phase function presents a reduction in memory requirement by ~20 times.

This study finally demonstrates that the appropriate use of approximated phase functions in a MC code can produce results reaching the one obtained if the exact phase functions were considered. The method M_{2} can be now applied in the case of inhomogeneous turbid media assuming that most of the drops considered are bigger than 10 /*μm* in diameter.

## 4. Application of the phase function approximation and comparison with experimental results for a hollow cone spray

#### 4.1. Monte Carlo simulation in a hollow cone spray

A typical planar Mie imaging experiment of a dilute, hollow cone spray generated by a Delavan pressure-swirl atomizer is simulated. The distribution of the extinction coefficient and of the dropsize come from experimental data obtained by the authors in [1] using laser induced fluorescence (LIF) and Phase Doppler Anemometry (PDA). The data is in the form of a 2D half plane and is shown in Fig. 11. The diameter of the droplets ranges from ~10 *μm* to ~75 *μm*.

The spray is assumed symmetrical and the full 3D structure is constructed in the model by rotating the data around the vertical axis. Each pixel on Fig. 11 represents a square area with 220 *μm* side, and the cubic cells in the MC model have the same side length. Figure 1 is a schematic of the MC simulation. The dimensions of the full simulated volume are 20 *mm* × 20 *mm* × 15 *mm*. The laser sheet (1 *mm* wide and 20 *mm* high) crosses the scattering medium in the middle of the spray. The wavelength of the laser is 532 nm. Drops are assumed spherical and non-absorbing with a refractive index of 1.4+0.0*i*. The detection area is one of the faces of the scattering volume parallel to the laser sheet (Fig. 1). The detector acceptance angle is *θ _{a}*=2.5°. With this angle a large number of photons are required to obtain good statistics. 5 billion photons are sent.

In each cell of the simulated volume, the average diameter of drops is given. Instead of taking the exact scattering phase function related to this average diameter, an approximate phase function is chosen using the method M_{2} described and verified previously.

#### 4.2. Results and discussion

Figure 12 shows the divergences between the experimental Mie image and the MC image. The laser light sheet enters on the left hand side of the image and leaves on the right. It can be seen that when all detected photons are included (Fig. 12(b)), the basic spray structure of the simulated image agrees well with the experimental image (Fig. 12(a)) even if some differences on the light intensity distribution can be noticed. These differences are explained by several factors: Firstly the restricted number of photons computed compromises the definition of the MC image. Secondly data used in the simulation are symmetrical around the spray axis, whereas real sprays of this type are known to be asymmetric by up to 15% in mass flow rate. Thirdly, MC data corresponding to the scattering coefficient have not been fully corrected from attenuation errors (described in introduction).

By taking into account all these factors the spatial resolution of MC images will reach the one obtain in experimental images allowing an accurate comparison of both images. Figure 12(c) is an image generated by numerically filtering Fig. 12(b) to include only singly scattered photons. It is seen that most of the detected photons are positioned on the left hand side of the spray (the side on which the laser sheet enters). Singly scattered light intensity is strongly reduced on the right edge of the spray image. However, as many of the scattering events are forwards scattering events with low angular deviations, the structure and direction of the light sheet is largely preserved, and the right part of the spray can be clearly observed in Fig. 1 a and b. It is seen that for single scattering the maximum number of detected photons per pixel is ~60 counts/pixel whereas with both single and multiple scattering the maximum is ~250 counts/pixel. Only 24% of the total number of detected photons has been singly scattered. The traditional assumption that all detected photons have been scattered only once, and carry information about single droplets only, is questionable. Multiple scattering occurring is dominant (76%) even for a spray assumed dilute and in which PDA measurements are possible. Multiple scattering occurs both along the light sheet and between the light sheet and the detector. Both cause undesirable errors in the detected signal if it is processed with the single scattering assumption. The magnitude of the error depends on the average deviation of the trajectory of the detected photons per scatter, and hence on the particle size distribution and the detector acceptance angle. Further work will concentrate on the quantification of the errors introduced by this multiple scattering.

## 5. Conclusion

A computationally efficient method to determine the scattering phase functions of different dropsizes has been presented and verified. This method allows saving consequent memory space (~20 times) and is of use in MC simulations of light propagation in complex inhomogeneous polydisperse sprays and aerosols in which both particle concentration and particle size distribution varies with location. The method is found to be valid for particle sizes from 10–200 *μm* diameter, with some differences with the rigorous method when droplets smaller than 10 *μm* are only considered. The method has been used to simulate a planar Mie imaging experiment in a hollow cone spray of moderate density. Simulated and experimental images have been compared and show that only 24% of the detected photons have been singly scattered. The use of the single scattering approximation is then questionable even for sprays assumed dilute. MC methods of type developed here are a necessary step in the development of new inverse techniques for improved optical diagnostics of such polydisperse inhomogeneous turbid media.

## Acknowledgments

The authors thank Thierry Mr Réveillé, Dr Thierry Girasole and Dr Claude Rozé for their help in the calculation of the average phase functions. This work was supported by EPSRC grant no. GR/R92653, by the Royal Society (project 15298) and NATO (project PST.CLG.979652).

## References and links

**
1
. **
E.
Berrocal
,
M.
Jermy
,
F.
Moukaideche
, and
I. V.
Meglinski
, “
Dense Spray Analysis using Optical Measurements and Monte Carlo simulation
,” presented at the 18
^{
th
}
Annual Conference on Liquid Atomization and Spray System-America, Irvine, CA, USA, 22–25 May
2005
.

**
2
. **
V.
Sicks
and
B.
Stojkovic
, “
Attenuation effects on imaging diagnostics of hollow-cone sprays
,”
Appl. Opt.
**
40
**
,
2435
–
2442
, (
2001
). [CrossRef]

**
3
. **
I. V.
Meglinski
,
V..P.
Romanov
,
D. Y.
Churmakov
,
E.
Berrocal
,
M. C.
Jermy
, and
D. A.
Greenhalgh
, “
Low and high orders light scattering in particulate media
,”
Laser Phys. Lett.
**
1
**
,
387
–
390
(
2004
). [CrossRef]

**
4
. **
R. P.
Meier
,
J. S.
Lee
, and
D. E.
Anderson
, “
Atmospheric scattering of middle UV radiation from an internal source
,”
Appl. Opt.
**
17
**
,
3216
–
3225
(
1978
). [CrossRef] [PubMed]

**
5
. **
T.
Girasole
,
C.
Roze
,
B.
Maheu
,
G.
Grehan
, and
J.
Menard
, “
Visibility distances in a foggy atmosphere: Comparisons between lighting installations by Monte Carlo simulation
,”
Int. Journal of Lighting Research and Technology
**
30
**
,
29
–
36
(
1998
). [CrossRef]

**
6
. **
R. F.
Bonner
,
R.
Nossal
,
S.
Havlin
, and
G. H.
Weiss
, “
Model for photon migration in turbid biological media
,”
J. Opt. Soc. Am. A
**
4
**
,
423
–
432
(
1987
). [CrossRef] [PubMed]

**
7
. **
I. R.
Abubakirov
and
A. A.
Gusev
, “
Estimation of scattering properties of lithosphere of Kamchatka based on Monte-Carlo simulation of record envelope of a near earthquake
,”
Phys. Earth Planet. Inter.
**
64
**
,
52
–
67
(
1990
). [CrossRef]

**
8
. **
M. C.
Jermy
and
A.
Allen
, “
Simulating the effects of multiple scattering on images of dense sprays and particle fields
,”
Appl. Opt.
**
41
**
,
4188
–
4196
(
2002
). [CrossRef] [PubMed]

**
9
. **
E.
Berrocal
,
D. Y.
Churmakov
,
V. P.
Romanov
,
M. C.
Jermy
, and
I. V.
Meglinski
, “
Crossed source/detector geometry for novel spray diagnostic: Monte Carlo and analytical results
,”
Appl. Opt.
**
44
**
,
2519
–
2529
(
2005
). [CrossRef] [PubMed]

**
10
. **
I. V.
Meglinsky
and
S. J.
Matcher
, “
Modeling the sampling volume for the skin blood oxygenation measurements
,”
Med. Biol. Eng. Comput.
**
39
**
,
44
–
50
(
2001
). [CrossRef] [PubMed]

**
11
. **
L.
Wang
,
S. L.
Jacques
, and
L.
Zheng
, “
MCML - Monte Carlo modelling of light transport in multi-layered tissues
,”
Computer Methods and Programs in Biomedicine
**
47
**
,
131
–
146
(
1995
). [CrossRef] [PubMed]

**
12
. **
I. V.
Meglinski
,
D. Y.
Churmakov
,
A. N.
Bashkatov
,
E. A.
Genina
, and
V. V.
Tuchin
, “
The Enhancement of Confocal Images of Tissues at Bulk Optical Immersion
,”
Laser Phys.
**
13
**
,
65
–
69
(
2003
).

**
13
. **
I. V.
Meglinski
,
V. L.
Kuzmin
,
D. Y.
Churmakov
, and
D. A.
Greenhalgh
, “
Monte Carlo Simulation of Coherent Effects in Multiple Scattering
,”
Proc. Roy. Soc. A
**
461
**
,
43
–
53
(
2005
). [CrossRef]

**
14
. **
B. T.
Wong
and
M. P.
Mengüç
, “
Comparison of Monte Carlo Techniques to Predict the Propagation of a Collimated Beam in Participating Media
,”
Numerical Heat Transfer
**
42
**
,
119
–
140
(
2002
). [CrossRef]

**
15
. **
H. C. van
de Hulst
, Light scattering by small particles (
Dover, N.Y.,
1981
).

**
16
. **
C.
Bohren
and
D.
Huffman
, Absorption and scattering of light by small particles (
Wiley, N.Y.,
1983
).

**
17
. **
L. G.
Henyey
and
J. L.
Greenstein
, “
Diffuse radiation in the galaxy
,”
Astrophys. J.
**
93
**
,
70
–
83
(
1941
). [CrossRef]