Image Grey Scale Based Volume Fraction Measurements of Solid Granular Particles

10TH INTERNATIONAL SYMPOSIUM ON PARTICLE IMAGE VELOCIMETRY - PIV13 Delft, The Netherlands, July 1-3, 2013

Image Grey Scale Based Volume Fraction Measurements of Solid Granular

Particles

Jari Kolehmainen1, Jouni Elfvengren1and Pentti Saarenrinne1

1_{Department of Energy and Process Engineering, Tampere University of Technology, Tampere, Finland}

jari.2.kolehmainen@tut.fi, jouni.elfvengren@tut.fi and pentti.saarenrinne@tut.fi

ABSTRACT

Fluidized beds are formed inside containers where the flow of fluid from the bottom of the container causes solid particles to be suspended in the fluid due to friction. As the terminal velocity increases the solid phase eventually gets fluidized and starts to behave in a more fluid-like manner. Typical industrial applications of fluidized beds range from energy production to the chemical industry. Fluidized bed combustion (FBC) has proved to be an efficient and low emission technology due to its intensely turbulent mixing of the solid and fluid phases. Further insights into this complex process are the topic of much currect research. Experimental measurements are required to validate the computational models, and for example, measurement integrated simulations. This work examines a relatively thin, almost 2 dimensional from the camera’s point of view, lab-scale fluidized bed partially filled with granular sand particles. While particle image velocimetry (PIV) can be used to measure the solid phase velocity, the volume fraction is hard to measure accurately. A novel volume fraction determination method for the solid phase is proposed in this work.

The proposed method is based on grey scale values obtained from images captured by a high-speed camera. These grey scale values are proportional to the incoming light intensity and can thus be regarded as a measure of the volume fraction. In the article [1] it was found that a logarithmic relation between the grey scale values and the solid volume fraction fits extremely well with the measured data. This work proposed a simple two-point calibration method for the brightest and the darkest average intensity values. Unfortunately this method only works if the bed is illuminated from behind. In dense suspensions, not enough light passes through to the front of the bed, which leads to poor PIV measurements. One solution for this problem is to illuminate the bed from the front, too. In consequence a more complex model has to be formulated for the relation between the grey scale values and the solid volume fraction. Furthermore, the calibration has to be performed for multiple solid volume fractions.

A five-parameter exponential model was developed to relate the greyscale values to the solid volume fraction. A set of calibration im-ages was used to fit the model parameters to the examined case. Grey scale histograms plotted from different calibration imim-ages showed that grey scale values follow a lognormal distribution. When the image brightness decreases, the measured grey scale distribution is found to get more centered around the average. The mean value of the cut lognormal intensity distribution was selected as a mapping parameter for the calibration curves. Based on the measurements, grey scale standard deviation was estimated to be roughly linearly dependent on the mean value of the cut lognormal distribution.

0 50 0 50 0.1 0.2 0.3 ε

(a) The real volume fraction field.

0 20 40 0 50 180 200 220 Intensity

(b) Image intensities of the real field.

0 50 0 500 0.2 0.4 ε

(c) Computed volume fraction field.

Figure 1: Method applied to an intensity field.

The proposed method was found to be applicable in the case of multiple light sources. The method can also cope with varying illumination conditions due to element-wise calibration that takes into account the changes in intensity within the image area. This is important since many laser-based light sources have a Gaussian intensity distribution, while beam shapers that would allow a uniform intensity distribution are expensive or limited by the power of the laser. In addition, the proposed method can take overexposure into account, producing a highly realistic model of the real situation.

The stability of the method was tested with a numerical simulation, as was its performance in static and dynamic cases. The simulations tests were focused on examining the effect of inferior parameters on the accuracy of the method. The method was found to work with

moderate errors in the calibration set. This is particularly important since exact calibration procedures may be time consuming or impossible to implement in practice. An example of a reconstructed volume fraction field can be seen in Figure 1.

The static case featured sand-filled tablets, in which the quantity of sand was known. The dynamic case consisted of a bubbling fluidized bed, where the mass balance was investigated. The results based on real data were promising, and due to its accuracy and flexibility the proposed method can be seen as a valuable aid in thin fluidized bed research.

1. Introduction

In a digital greyscale image, the light intensity values are stored as discrete counts. The commonly used 8-bit grayscale consists of 256 shades of grey for different intensity values. As a standard, the minimum value 0 equals a totally black pixel and the maximum value 255, is used fora totally white, overexposed pixel. The greyscale values are usually generated in a CCD or CMOS camera. In a single pixel of the camera sensor, the photons of incoming light cause the electrons of the semiconductor material to rise from the valence to the conduction band. As a result the number of electrons observed as a charge is linearly proportional to the intensity of the light. However, most of the basic cameras utilise a brightness gradation called gamma correction for improved visualization of the image, so the relationship between a greyscale value and the actual sensor charge is not necessarily linear. In professional cameras the linear intensity to greyscale data is usually available and is often stored as a default [4]. For example, in the Photron Fastcam SA5, the brightness is not converted as long as none of the specific look-up table operations have been selected. The above-described linear behavior of the greyscale values forms the basis of the volume fraction measurements, although distortion caused by overexposure also has to be taken into account.

One might assume that imaging-based volume fraction measurements from a thin fluidized bed would require very homogenous il-lumination conditions. However, in practice a uniform light intensity distribution for the whole imaging area is often difficult, and certainly very expensive to achieve. This is because many light sources tend to concentrate the light in a somewhat Gaussian intensity distribution. It would be possible to utilize equalizing optics to obtain a uniform intensity distribution, but such devices are expensive. An alternative elemental approach, which allows uneven illumination conditions, is presented here.

2. Measurement setup

The experiments were done on a small thin fluidized bed. The bed’s walls were made of transparent plastic, and the ”‘sand”’ used in the bed consisted of small spherical glass particles with an average diameter of 0.2mm. The high speed camera was positioned facing the bed. The light was supplied by two pulsed Cavitar HF diode lasers with wavelengths of 810nm, one in front of the bed and one behind it relative to the camera. The laser light intensities were controlled by limiting the pulse width. In addition, a bandpass filter was used in order to decrease the noise in the pictures. The filter was designed for an 810nm wavelength and had a bandwidth of 10nm. Figure 2 shows a schematic diagram of the measurement setup.

Optical fibre Optical fibre Fluidized bed Diffuser Bandpass filter Camera

Figure 2: Schematic diagram of the measurement setup.

Before the actual measurements can be taken from the fluidized bed, the illumination setup has to be calibrated using calibration tablets. These have a known thickness which is directly related to the volume fraction in the 2D case. In this work, eight distinct tablets were used for calibration. Multiple images were obtained from each tablet, since shaking or moving the tablet between the snapshots can change the mean intensity slightly.

3. Model

3.1 Conditional Image Intensity

The image intensity in an image of a fluidized bed varies within the bed depending on the locations of the particles and the illumination conditions. In this work, the main interest is not the exact particle locations, but rather the mean amount of particles in a given location.Therefore, it is practical to treat image intensity as a random variable. The image intensity ranges from 0 to 255 depending on the voltage of the sensor in the camera. Light intensities larger than the maximum observable light intensity were treated as 255 and those which were lower than the minimum observable light intensity were treated as 0.

Experiments revealed that the image intensity behaves like a lognormal distribution [2] in a given snapshot of a fluidized bed when image intensities are not near the boundaries of the camera’s observable intensity range. The parameters of the distribution depended on the illumination conditions and the volume fraction of particles in the given location. However, when the intensity values approached the maximum value of the observable intensity range there was a large spike at the maximum intensity value. This phenomenon was

not observed at the minimum intensity value. Because of the partial overexposure of images, a cut lognormal distribution was used as a model for the image intensity with a known volume fraction and uniform illumination conditions.

p(I|µ, σ) = U (0, IM)φ(I|µ, σ) + Φ(µ − IM|µ, σ)δ(I − IM), (1)

where δ is Dirac’s delta function, µ and σ the parameters of normal distribution, U (0, IM) the constant distribution from 0 to IM, φ is

the density function of the lognormal distribution and Φ the cumulative distribution function of lognormal distribution. The parameter IMstands for the upper limit, which is 255 in this study.

Given an image with image intensities Ii, the parameters µ and σ can be found by expectation maximization:

(µ∗, σ∗) = argmax_µ,σp(µ, σ|I) (2) = argmax_µ,σ

_∑

_∑

_∑

Problem (2) can be solved by conventional nonlinear optimization algorithms such as SQP. A good initial guess for the problem can be obtained by fitting a lognormal distribution to the data.

The mean and variance of a lognormal random variable with parameters µ and σ are given by

E(I|ε) = exp(µ + σ2/2) (5)

var(I|ε) = (exp(σ2) − 1)E(I|ε)2. (6)

Moreover, it was experimentally found that the mean intensity (5) is related to the volume fraction by the exponential model

E(I|ε) = a1+ a2exp(a3ε) + a4exp(a5ε), (7)

which takes the influence of simultaneous front and back illumination into account. Parameters aiare experimental and sensitive to the

illumination conditions.

Interestingly, the parameter σ was not too sensitive to changes in illumination conditions or volume fraction, which in the light of Equation (6) suggests that the standard deviation of image intensity is proportional to the mean image intensity. Hence, it holds that

var(I|ε) = γ2E(I|ε)2, (8)

where, based on calibration data, gamma was, approximated as a constant γ ≈ 0.19.

Combining Equations (5), (6) and (7) gives a relation between the distribution parameters and mean intensity, which can be computed from volume fraction by Equation (7). However, the sample mean ˆIcomputed from the picture is not the same as the mean in Equation (5), but rather the mean value of the cut lognormal distribution. The sample mean is connected to the value (5) by

ˆ

I− E(I) = IM(Φ(IM|µ, σ) − 1) −

Z∞

IM φ(I|µ, σ)IdI.

(9)

Equation (9) can be solved by, for instance , fixed point iteration, since µ and σ are known functions of E(I). The calibration parameters of model (7) were fitted to the mean values E(I) solved from Equation (9).

3.2 Calibration procedure

First, a calibration image with a known volume fraction is divided into elements of size dx, dy in pixels. The element size is typically selected to equal the size of the PIV interrogation window. The element’s mean intensity is computed from the greyscale values of the pixels inside each element using Equation (5). The distribution parameters are found by solving the optimization problem (2). This is repeated to a set of calibration images including different volume fractions.

Next the parameters of model (7) are computed for each element using the Levenberg-Marquard method. These parameters are later used in the volume fraction computations of the fluidized bed model. The intensity deviation inside the imaging area is taken into account by the elemental approach of the calibration procedure. It should be noted, though, that in very small element sizes the degree of uncertainty increases due to the smaller sample size of the pixels inside the element. To counteract this effect, several images are taken from a single calibration tablet.

After the calibration, the fluidized bed model is placed in the exact location where the calibration images are taken. Camera options (calibration, aperture, focus, resolution etc.) and illumination conditions (intensity, distances) must remain unchanged during the calibration and the actual recordings. Furthermore, solid optical properties, such as reflections from the surfaces of the calibration tablets and the model, should remain quite similar. These conditions ensure that the parameters obtained from the calibration are applicable in the volume fraction calculation of the actual model.

3.3 Bayesian Solution

Even though the relation between image intensity and volume fraction is simple, the inverse problem becomes ill-posed when the volume fraction gets large. Moreover, the variance of the mean image intensity of an element increases as the size of the element decreases. This can cause sudden jumps in the maximum image intensity and unwanted and unrealistic oscillations in the solution when simple inversion methods, such as Newton iteration are used.

One way to cope with the above mentioned behavior is to use Bayesian statistics. In particular, this approach enables one to use the variance data to its full extent. The posterior density of the volume fraction is given by Bayes’ rule

πi(ε|I) ∝ p(ε)pi(I|ε), (10)

where p(ε) is a prior density and i refers to the corresponding element. In this work, a noninformative prior was used, which works well on large to mid-sized elements. Small elements might need regularization, which could be implemented in to the prior distribution. The chosen prior was

p(ε) = U (ε|εm, εM), (11)

where εmis the lower packing limit and εMis the upper packing limit.

In this work, the chosen estimator of the measured volume fraction was the conditional expected value of the posterior density. The conditional expected value is given by

˜εi= E (ε|I) = Z_εM εm pi(I|ε)dε −1Z_εM εm pi(I|ε)εdε. (12)

The variance of the estimator ˜ε was estimated by

var (˜εi) = ZεM εm pi(I|ε)dε −1ZεM εm pi(I|ε) (ε − ˜ε)2dε. (13) 4. Testing of method 4.1 Simulation tests

The stability of the method was tested by altering the calibration set. This alteration was carried out by adding normal noise with known mean and variance to each intensity value in the calibration set. This introduced erroneous calibration parameters which were used to compute a reconstruction of a user specified volume fraction field. The computed field was compared with the exact volume fraction field, and the expected value of the RMSE (Root Mean Square Error) was computed using a Monte Carlo-simulation. The sample size in the Monte Carlo-simulation was 5000 samples for each error condition. This shed some light on the method’s behavior when the calibration process is not exact. Results are shown in the Figure 3.

These results were promising as the RMSE increased in a linear manner as the mean error increased. Moreover, the variance of the noise only seemed to have a minor effect on the accuracy of the method, which is a highly desirable feature.

4.2 Test with a bubbling fluidized bed

In a bubbling fluidized bed the grate pressure is too low to cause fluidization and, air escapes from the grate by forming large bubbles. An example of time series taken from a bubbling fluidized bed is shown in Figure 4. The bed used in this experiment was 6.5mm thick and 100mm wide. These dimensions forced the bubbling bed to behave in a two dimensional manner. The diode lasers were positioned behind and in front of the bed. The front laser had a pulse width of 0.2µs while the more dominant back laser had a pulse width of 0.59µs. Calibration was done using eight calibration tablets with thicknesses of 0.75mm, 1.00mm, 1.25mm, 1.50mm, 2.00mm, 3.00mm, 4.00mm and 6.00mm. The frame rate was set to 500fps with full exposure and the pinhole ratio was 2.8. The initial particle volume was measured as 40mL ± 0.1mL.

The proposed lognormal based method was compared to logarithmic model based on direct use of Beer’s law and a Bayesian method similar to one proposed in this study with exception that the lognormal distribution was replaced by a normal distribution. The lognor-mal model is described in detail in the thesis [3].

The results of the experiment are shown in Figure 5a. As one can see from the figure, the proposed Bayesian approach has better accuracy than the logarithmic model in each frame. Moreover, the lognormal distribution was clearly better than the normal distribution based model in this test.

0 0.5 1 1.5 2 0 1 2 3 4 0 0.2 0.4 0.6 0.8 1 1.2 1.4 ε σ RMSE

Figure 3: RMSE of the volume fraction field with different error conditions in the calibration set.

4.3 Test with a turbulent fluidized bed

In a turbulent fluidized bed the grate pressure is larger than in a bubbling bed, which results in a continuous fluidization. An example time series of a turbulent bed is shown in Figure 6. The experimental setup was the same as with the bubbling bed described in the previous section. In the turbulent case, the exact volume at the image domain is unknown. However, one can investigate how well the different methods satisfy mass conservation. Fluxes were evaluated from the corresponding volume fraction fields coupled with PIV measurements. The relative error of mass conservation was computed by

ε = kR V∂ε∂t+ ∇ · (εu) dV k R VεdV , (14)

where integration is done over the whole image domain. Moreover, the discretization of Equation (14) is done with respect to the PIV interrogation scheme.

The results of the comparison are shown in Figure 5b. The proposed lognormal distribution based method did have larger error than the normal-distribution-based model. However, the test setup only gives an insight into the actual error, which is not the same as the relative error. Moreover, the PIV fields are not error-free, and hence may cause additional bias in the actual order of the methods. Thus the results of the turbulent bed case are not accurate and only give some insight into the errors of the different methods.

(a) Frame at 0ms (b) Frame at 10ms (c) Frame at 20ms

(d) Frame at 30ms (e) Frame at 40ms (f) Frame at 50ms

Figure 4: Time series taken from a bubbling fluidized bed between 10ms.

0 1000 2000 3000 4000 5000 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Frame Relative error

(a) Relative error

0 100 200 300 400 500 −4 −2 0 2 4 6 8 10 Frame log(ε)

(b) Logarithm of relative error

Figure 5: Method applied to different experiments. Logarithmic model is shown with red curve, normal distribution based model with

the green curve, and lognormal distribution based model with blue curve.

5. Conlusions

The advantage of using multiple lasers rather than just one, is that the quality of the measured PIV field increases greatly in the dark parts of the fluidized bed due to the introduction of the second laser. Moreover, the omission of this effect in the modelling resulted in increased error as can be seen from the comparisons.

(a) Frame at 0ms (b) Frame at 10ms (c) Frame at 20ms

(d) Frame at 30ms (e) Frame at 40ms (f) Frame at 50ms

Figure 6: Time series taken from a turbulent fluidized bed between 10ms.

The proposed method performed best in the case of the bubbling fluidized bed. The test suggested that the additional computational cost of using lognormal distribution instead of normal distribution is justified. Moreover, the Bayesian approach performed distinctly better than the logarithmic model based directly on Beer’s law. The stability of the proposed methods was investigated using a Monte Carlo-simulation, which suggested that the calibration can tolerate moderate errors and even biased calibration tablets, with linear increase in the error in the solution.

From the computational point of view, the Bayesian approach is slower than the logarithmic model, but it parallelises quite well. This suggests that the method could be implemented in a GPU quite efficiently, which would enable processing of large data sets.

REFERENCES

[1] G. Grasa and Abanades J.C. A calibration procedure to obtain solid concentrations from digital images of bulk powders. Powder Technology, 114(1–3):125–128, January 2001.

[2] Norman L. Johnson, Samuel Kotz, and N. Balakrishnan. Continuous Univariate Distributions, volume 1. Wiley and Sons, 2nd edition, 1994.

[3] J. Peltola. Dynamics in a circulating fluidized bed: Experimental and numerical study. Master’s thesis, Tampere University of Technology, June 2009.