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:

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(ξ)/μ_ext with μ_ext =N ^* σ_ext. N is the number density of scattering particles and σ_ext 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=μ_sca+μ_abs where μ_sca is the scattering coefficient and μ_abs is the absorption coefficient) and the albedo equal 1 (Λ = μ_sca/(μ_sca + μ_abs)). 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 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 θ_s defined between 0 and π is found from the inverse CPDF of f by θ_s = CPDF^-1(ξ) (where ξ is a random number between 0 and 1). The azimuthal scattering angle φ_s is uniformly distributed between 0 and 2π. 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.

Table 1. Scattering turbid media types and required optical properties

View Table

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 μ_ext is in this second case calculated from the average extinction cross section σ¯_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.

 

Fig. 1. Illustration of Monte Carlo modeling in inhomogeneous polydisperse scattering medium

Download Full Size | PDF

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.

 

Fig. 2. Normalized Log-Normal distributions of particle size for different d¯ and σ.

(a): d¯ = 40 μm and σ = 32 μm (80% of d¯), (b): d¯ = 40 μm and σ = 4 μm (10% of d¯),

(c):d¯ = 5 μm and σ = 4 μm (80%ofd¯), (d): d¯ = 5 μm and σ = 0.5 μm (10%ofd¯).

Download Full Size | PDF

Two methods for representing the local scattering phase function in polydisperse homogeneous media were tested. The method one (M₁) 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₁. In the case of infinitesimal size bin width of the dropsize distribution M₁ 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₂), 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 ₁) of a drop of diameter d ₁ 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₂ would be equal to method M₁. The accuracy of M₂ 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.

 

Fig. 3. Representation of 25 Cumulative Probability Density Functions (CPDF) calculated from Mie theory. Each CPDF represents a class of drop size. Scattering particle diameters range from 2 to 200 μm.

Download Full Size | PDF

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).

 

Fig. 4. The scattering medium is a single homogeneous cube of L = 50 mm. The source S is a cylindrical laser beam.

Download Full Size | PDF

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.0i). The droplets are spherical and non-absorbing with refractive index 1.4+0.0i and the scattering CPDF are calculated from the Mie theory. Scattering and extinction coefficients are then equal and the simulations are run with μ_sca fixed to 0.12 and 0.24 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₂ are compared to the results obtained from M₁ 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₁ 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.12mm ^-1 and 0.24 mm ^-1 giving respectively an average number of scatters per photon of ~6 and ~12 scatters. Figure 5(a) (μ_sca =0.12 mm ^-1) and (b) (μ_sca =0.24 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 μ_sca 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.

 

Fig. 5. Intensity distribution for forward (front face) and backward (back face) light scattering at different scattering coefficient. μ_sca =0.12 mm ^-1 for (a) (c) (e) (g) and μ_sca =0.24 mm ^-1 for (b) (d) (f) (h). The detection acceptance angle θ_a = 90°. 100 million photons are sent. The intensity scale represents the number of photons detected per pixel. Pixels are 25 μm sides. [Media 1] [Media 2] [Media 3] [Media 4]

Download Full Size | PDF

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 μ_sca = 0.24 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 θ_a = 90°). As found in other simulations [9], the acceptance angle can be used to optimize the ratio of singly to multiply scattered photons detected.

 

Fig. 6. Intensity distribution for forward light scattering (front face). μ_sca =0.12 mm ^-1 for (a) (c) and μ_sca =0.24 mm ^-1 for (b) (d). The detection acceptance angle θ_a = 5°. 100 millions of photons are sent. The intensity scale represents the number of photon detected per pixel. Each pixel is square with 25 μm sides. [Media 5] [Media 6]

Download Full Size | PDF

In Fig. 6 the detection acceptance angle θ_a is reduced to 5°. μ_sca=0.12 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 θ_a = 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 θ_s less than 5° and shows the high intensity of scattering in the forward direction in Mie scattering processes.

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₂) 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₁. Corresponding images have been generated by applying M₂. 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₁ to the image generated with M₂. This ratio is plotted (Fig.7) for μ_sca= 0.12 mm ^-1 and μ_sca= 0.24 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 θ_a = 90°. It is seen from Fig.7 that the ratio M₁/M₂ (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.

 

Fig. 7. Intensity ratio, along a beam profile, between the 2 MC methods. Forward scattering is in (a) and (b) and backscattering in (c) and (d). Both single and multiple scattering are detected together. Detector acceptance angle θ_a = 90°.

μ_sca =0.12 mm ^-1 for (a) (c) and μ_sca =0.24 mm ^-1 for (b) (d).

Download Full Size | PDF

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₁/M₂ 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₂ match the results from M₁ 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₁/M₂ intensity ratio noisier for μ_sca= 0.24 mm ^-1 Fig. 8(b) than for μ_sca = 0.12 mm ^-1 Fig. 8(a).

 

Fig. 8. Intensity profile ratio between the 2 MC methods for forward light scattering. Single and multiple scattering are detected with an acceptance angle θ_a = 5°. μ_sca =0.12 mm ^-1 for (a) (c) and μ_sca =0.24 mm ^-1 for (b) (d).

Download Full Size | PDF

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 θ_a=90° (Fig. 9(a)), few photons are detected giving strong statistical fluctuations in the resulting ratio M₁/M₂. If however θ_a=5° (Fig. 9(b)), no singly scattered photons are detected for either method and a flat line is plotted for X > 40 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.

 

Fig. 9. Intensity profile ratio between the 2 MC methods for forward light scattering with μ_sca =0.12 mm ^-1. Single scattering only is detected with an acceptance angle θ_a = 90°(a) and θ_a = 5°(b). No singly scattered photons are detected with either method at X>40mm in (b).

Download Full Size | PDF

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₁ than M₂ 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₁ than in M₂.

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₂ are in excellent agreement with the rigorous method M₁, for droplet sizes over the range 2–200 μm of particle size, for different θ_a=5° and 90° and with μ_sca varying over a factor of 2.

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.

 

Fig. 10. Intensity profile ratio between the 2 MC methods for μ_sca =0.12 mm ^-1. Single and multiple light scattering are detected on the front face with an acceptance angle θ_a = 90°. The log-normal distributions of particle size are characterized by: (a): d¯ = 40 μm with à = 4 μm (10%), (b): d¯ = 5 μm with σ = 4 μm (80%), (c): d¯ = 5 μm with σ = 0.5 μm (10%).

Download Full Size | PDF

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₁/M₂ 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₁/M₂ < 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₂ 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₂ (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₂ 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.

 

Fig. 11. Extinction coefficient (a) and droplet diameter (b) through the central plane of a hollow cone spray. By rotating the data around the central vertical axis MC input data are generated in 3D.

Download Full Size | PDF

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.0i. 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₂ 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).

 

Fig. 12. Comparison between the experimental (a) and MC (b,c) images for a planar Mie imaging of a hollow cone spray. (c) is generated from the singly scattered detected photons only.

Download Full Size | PDF

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.