Riley Re 2009 Thesis

Image analysis of turbulent mixing after a fractal Cantor grid Richard E Riley May 2009

Supervisor: Dr F. Nicolleau Thesis submitted to the University of Sheffield in partial fulfilment of the requirements for the degree of Master of Engineering

Abstract The turbulent mixing of a marker species in a flow of water is investigated using image analysis. A transient mixing regime is established comprising a homogeneous upstream concentration, interrupted by a grid. The development of the downstream concentration field is recorded in a sequence of digital images. Grid pairs are tested, each consisting of a fractal grid in the form of extruded Cantor dust (generations one to four), and a non-fractal control grid of equivalent flow area. The fractal dimension of the species interface is calculated over the sequence using the boxcounting method. Statistical algorithms and a method of wavelet analysis are developed and applied as measurements of mixing. A novel method of species identification using pixel saturation is employed and is found to produce clearer results, independant of shadow effects and non-uniform lighting. Flow through higher generation grids is observed to be near-laminar. Where the flow is turbulent, no special behaviour is observed in terms of the fractal dimensions of the species interface. Wavelet analysis alludes to increased energy at the smaller scales of the fractal grid concentration fields. Improvements to the method are discussed and an experiment is proposed involving a stationary mixing regime for more appropriate analysis. R´ esum´ e On ´etudie le m´elange turbulent d’une teinture dans un courant d’eau par analyse d’images. Le r´egime transitoire du m´elange se compose d’une concentration homog`ene en amont d’une grille. Le d´eveloppement du champ de concentration en aval de la grille est enregistr´e dans une s´equence d’images num´eriques. L’´etude concerne des grilles par paire; une grille fractale en forme extrud´ee de la poussi`ere de Cantor, et une grille de contrˆ ole avec une aire d’´ecoulement ´equivalente. On calcule la dimension fractale de l’interface d’esp`ece `a l’aide de la m´ethode de boxcounting. On d´eveloppe des algorithmes statistiques et une m´ethode d’analyse d’ondelette que l’on applique comme mesure de m´elange. On utilise une nouvelle m´ethode d’identification, bas´ee sur la saturation de pixel, qui fournit des r´esultats

i

plus clairs et ind´ependants des effets d’ombres o` u d’´eclairage non-uniforme. On observe l’´ecoulement quasi-laminaire `a travers les grilles de g´en´erations sup´erieures. Dans le cas o` u l’´ecoulement est turbulent, on n’observe pas de phenom`ene sp´ecial en termes de dimension fractale d’interface. L’analyse d’ondelette fait allusion `a une augmentation de l’´energie aux ´echelles inf´erieures du champ de concentration des grilles fractales. On propose des am´eliorations `a la m´ethode ainsi qu’une exp´erience qui comprend un r´egime stationnaire de m´elange pour l’analyse plus comp´etent.

ii

C(x, y)

Marker species concentration field

-

CoV

Coefficient of variance

-

D

Dimension

-

E

Exposure

-

G(λ)

Global wavelet energy as a function of scale

-

(m, n)

Sampling point coordinates in discretised space

-

N

Number of covering objects

-

r

Eddy scale diameter

Re

Reynolds number

m -

S(m, n) Pixel saturation map

-

(x, y)

Rectangular coordinates in continuous 2-D space

-

S0

Threshold saturation

-

Np

Number of sampling points

-

λ

Scale/wavelet scale length

m

iii

Contents 1 Introduction

1

2 Theory

3

2.1

2.2

2.3

2.4

The structure of turbulence . . . . . . . . . . . . . . . . . . . . . . . . .

3

2.1.1

The energy cascade . . . . . . . . . . . . . . . . . . . . . . . . . .

4

2.1.2

Practical applications . . . . . . . . . . . . . . . . . . . . . . . . .

6

Mixing levels and mechanisms . . . . . . . . . . . . . . . . . . . . . . . .

6

2.2.1

Measuring mixing . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

Wavelets . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

9

2.3.1

The wavelet transform . . . . . . . . . . . . . . . . . . . . . . . .

9

2.3.2

The global wavelet energy spectrum . . . . . . . . . . . . . . . . .

11

Fractals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

12

2.4.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

12

2.4.2

Fractal dimension . . . . . . . . . . . . . . . . . . . . . . . . . . .

13

2.4.3

Box-counting . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

16

3 Procedure 3.1

3.2

17

Experimental procedure . . . . . . . . . . . . . . . . . . . . . . . . . . .

17

3.1.1

Apparatus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

17

3.1.2

Flow conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . .

17

3.1.3

Image capture . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

19

3.1.4

Dye injection . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

20

3.1.5

The grids . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

22

Numerical analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

22

3.2.1

Box-counting . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

26

3.2.2

Wavelet analysis . . . . . . . . . . . . . . . . . . . . . . . . . . .

26

3.2.3

Statistical analysis . . . . . . . . . . . . . . . . . . . . . . . . . .

26

iv

4 Results

27

4.1

Statistical properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

28

4.2

Global wavelet energy spectrum . . . . . . . . . . . . . . . . . . . . . . .

29

4.3

Fractal dimension . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

30

5 Discussion

32

5.1

Four stages of development . . . . . . . . . . . . . . . . . . . . . . . . . .

32

5.2

Statistical properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

34

5.3

Wavelet analysis - the scales of segregation . . . . . . . . . . . . . . . . .

35

5.4

Fractal dimension . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

36

5.5

Grid comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

37

5.6

Method review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

38

6 Conclusion and Recommendations

41

References

44

Appendices

46

A Image sequences

46

B Ancilliary information

55

C Algorithms

57

C.1 Time meta-data extraction . . . . . . . . . . . . . . . . . . . . . . . . . .

57

C.2 Image conversion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

60

C.3 Box-counting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

62

C.4 Box-counting batch code . . . . . . . . . . . . . . . . . . . . . . . . . . .

68

C.5 Statistical analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

70

C.6 Global wavelet energy spectrum . . . . . . . . . . . . . . . . . . . . . . .

73

v

Acknowledgements The author would like to thank the following individuals: Dr Franck Nicolleau and Dr Andrzej Nowakowski for their continued input and enthusiasm throughout the project. Mr Wayne Oxley for his invaluable contribution to the construction of the apparatus. Dr Mike Green and his colleagues at the Low Carbon Combustion Centre for the extended loan of the digital camera. Mr Malcolm Nettleship for the supply of a flow meter.

vi

1

Introduction

Turbulence is an extremely important phenomenon of fluid dynamics which plays a role in many industrial processes. The topology of turbulence is of particular interest; the eddies which characterise turbulent flow exist across a range of physical scales pertaining to the transfer of energy from decaying larger scales to be dissipated viscously at the smallest scales. The theory of fractal geometry also describes self-similar repeating scales. The fractal nature of turbulence was first examined in detail by Sreenivasan and Meneveau in 1986 [7]. Much subsequent work has been conducted by Catrakis and Bond, most significantly, for this investigation at least, with regards to the numerical box-counting method as a tool for calculating the fractal dimension of the turbulent/non-turbulent interface. Turbulence is highly dissipative and plays a large role in fluid mixing, which has a number of important applications such as combustion, chemical synthesis, petro-chemical production and water & sanitation, as well as environmental problems such as contaminant transfer in ground water [8]. This investigation examines the development of turbulence, and its effect on mixing, by observation of a marker species in a flow of water through a grid. Analysis of the resulting concentration field is undertaken using digital image analysis. Box-counting, wavelet analysis and statistical methods are employed to examine mixing and scale development in the flow. Mimicking the fractal nature of turbulence, fractal Cantor grids are compared to nonfractal control grids of equivalent flow area. The investigation aims to examine the effect of the Cantor grids on mixing and the development of turbulence. Particular attention is paid to the fractal dimension of the dye interface, and the statistical properties and development of scales in the concentration field. The objectives of the investigation are as follows: ˆ To progress towards a thorough laboratory procedure which is fully controlled.

1

ˆ Develop methods of analysis which are standardised and comparable across all

results. ˆ Investigate the effects of fractal and non-fractal grids on mixing and scale develop-

ment. ˆ Make recommendations for future work.

2

2 2.1

Theory The structure of turbulence

A turbulent flow is defined by the dominance of inertial forces over viscous forces which would otherwise tend to damp the response to perturbations in the flow field. The result is instability of the flow and a sensitivity to small alterations of the boundary conditions. Turbulence is characterised by a number of important properties. Randomness in time and space Turbulence is a random process, comprising a large number of degrees of freedom. The result is unpredictable on a detailed level; however, the statistical properties are repeatable, allowing analysis on the macro scale [9]. A wide range of scales A turbulent flow may be viewed as a composition of eddies of different sizes, or scales [12]. Energy cascades from the large driving scales to the smallest scales, where it is dissipated viscously. These scales are superimposed in time and space, accounting for the complexity of a turbulent field. A brief look at the non-linear term of the Navier-Stokes equation indicates the origin of these scales in the mathematical model. ∂ 1 ∂ ∂ ∂2 Ui + Ui Uj = − P +ν Ui ∂t ∂xj ρ ∂xi ∂xi ∂xj | {z }

(2.1)

non-linear

Highly dissipative Turbulence is highly dissipative in terms of heat as well as energy, and dispersive in terms of mass transfer, which is beneficial to fluid mixing [9]. Three-dimensional Turbulence is the creation of vortices, which interact in 3-D space [9].

3

Figure 2.1: An illustration of the energy cascade. 2.1.1

The energy cascade

The flow of energy from large to small scales is known as the energy cascade. Kinetic energy is introduced at the driving scales which, due to instability, decay into smaller sub-scales. These sub-scales decay into yet smaller scales as illustrated in figure 2.1. This process continues until viscous forces become dominant and energy can be dissipated viscously. This theory was proposed by Richardson (1922) and it would be churlish not to include his very effective poetic summary of the process: Big whorls have little whorls, which feed on their velocity; and little whorls have lesser whorls, and so on to viscosity. The smallest scale present in a turbulent flow field - that at which energy is dissipated viscously - is known as the Kolmogorov length scale, η, after Andrei Kolmogorov (1941). Kolmogorov posed that this scale corresponds to a local Reynolds number Reη ≈ 1, based on the characteristic dimension of the eddy. This signifies the point at which viscous forces become locally dominant. It is possible to relate the Kolmogorov scale to the driving scale, L, by considering the transfer of energy through the cascade. 4

The driving energy of the flow is of the order of the square of its mean velocity. E = O < u2 >

(2.2)

A characteristic time, T , can be defined as follows.

T =√

L < u2 >

(2.3)

Therefore, the energy dissipated over time T may be defined as the energy dissipation rate, ²¯. 3

< u2 > < u2 > 2 ²¯ = = T L

(2.4)

At a given scale, r, the energy flux, ²¯, may be defined.

²¯ =

u3r r

(2.5)

By rearrangement, the velocity at r may be expressed as:

ur =

√

<

u2

>

h r i 13 L

(2.6)

When Reη ≈ 1 it can be shown that: h η i 43 L

=

ν √ L < u2 >

(2.7)

It is therefore possible to estimate the Kolmogorov length scale with some knowledge of the driving scales. In general, this knowledge is obtained from direct velocity measurements.

5

2.1.2

Practical applications

The dispersive effect of turbulence has implications in mixing. Flow mixing is important in a number of engineering applications as well as natural processes. Many processes require the mixing of two or more species. For example, combustion relies on the mixing of fuel and air, while the production of a soft drink requires the mixing of the ingredients. Understanding and improving mixing can improve the efficiency of these processes. This investigation concerns a continuous mixing problem. Continuous mixing offers some obvious advantages over batch processes. Continuous mixing systems require less space, are more easily controlled for improved consistency, and require reduced labour [11].

2.2

Mixing levels and mechanisms

In the most basis case, mixing is passive and does not affect the flow dynamics, which may be considered separately. Dimotakis [2] labelled this as level 1 mixing. Level 2 mixing describes a regime where the flow field properties are not independent of the mixing process. Such mixing may occur between fluids of a different density, temperature or salinity. Level 3 mixing occurs when a change in the fluid results from chemical reaction. A typical example is combustion. Three basic processes of fluid mixing exist: stretching & folding, diffusion and breakup [10]. The turbulent mixing of species can be broken down into 3 mechanisms, as described by Eckart in 1948; entrainment, kinematic stirring (dispersion) and molecular mixing. Entrainment At the largest scales of turbulence, fluid in non-turbulent regions is entrained into the turbulent flow field.

6

Kinematic stirring Intermediate scales are responsible for stirring (folding and stretching), thus creating a large species interface and accelerating the mixing process. Molecular mixing The final stage of mixing occurs at the small scales and is associated with viscous dissipation and molecular diffusion [1]. Molecular and viscous diffusion rely on a concentration gradient. Fick’s first law states that the mass flux per unit area associated with molecular diffusion is proportional to the concentration gradient. −−→ − → JA = −δAB · gradCA

(2.8)

− → The mass flux vector of fluid A per unit area is denoted JA , δAB is the coefficient of diffusion of species A in fluid B, and CA is the concentration of species A. In level one mixing, viscous diffusion is analogous to molecular diffusion, also relying on a gradient of concentration [2]. Mass transfer occurs from regions of high concentration to regions of low concentration, as indicated by the negative sign. 2.2.1

Measuring mixing

One may easily define a mixed and unmixed state, but intermediate states are more problematic. Kukukova et al. [4] have proposed three dimensions to describe mixing. During this investigation, the mixing of one marker species is considered. Concentration of the marker species is defined as C. Intensity of segregation The extent to which a species can be considered separate is described by the intensity of segregation. A concentration field with an intensity of segregation of unity contains only distinct homogeneous regions of a species. A value of zero indicates that the entire field contains a homogeneous concentration. The intensity of segregation can be defined as the coefficient of variance, v u ¶2 Np µ u 1 X Ci − Cmean t CoV = Np i=1 Cmean 7

(2.9)

where Np is the number of sample points and Ci is the marker concentration at point i. The coefficient of variance does not distinguish between a total absence of a species and a homogeneous mix of a species. Exposure Exposure may be considered a description of the length and sharpness of the species interface. According to Kukukova et al. (2009) it is “analogous to a simplified calculation of the rate of mass transfer across the interface ”[4], as implied by Fick’s law. A simplified and normalised calculation of exposure, E, is used in this investigation, Np Na 1 XX 1 ∼ E= |Ci − Cj | Np i=1 j=1 2

(2.10)

where Na is the number of adjacent sample points (2, 3 or 4 for a square sampling regime) to the current sampling point i. Scales of segregation Break-up can be described by the scales of segregation. The scales of segregation are difficult to define for a large quasi-continuous field such as a digital image. Wavelet analysis provides a useful tool and is discussed in more detail in section 2.3. A priori, in a state of homogeneous concentration (or a ‘perfect’mix), all three dimensions tend to zero. During this investigation, the concentration field, C(x, y), is represented by a quasicontinuous saturation map, S(m, n), calculated from the digital images (see 3.2) - all analysis is applied to this digital representation.

8

2.3

Wavelets

Wavelet analysis allows the decomposition of signal data into space and local scale information, which is useful in the analysis of scales and provides a measure of the scales of segregation in a concentration field. The one-dimensional wavelet transform is introduced and the global wavelet energy spectrum is defined as a method of consolidating the scale information into an easily comparable form. 2.3.1

The wavelet transform

The wavelet transform decomposes a function into a combination of wavelet packets which are translated and scaled versions of a carefully chosen mother wavelet. The wavelet transform is closely related to the Fourier transform, which decomposes a signal into trigonometric functions. The translation and scale of the packets provides information on a function’s scale composition in local space. The mother wavelet is an oscillating function of finite length, subject to certain conditions [3]. To be oscillatory, the mother wavelet must meet the following condition, Z

∞

ψ(y)dy = 0

(2.11)

−∞

which implies that it is negative as often as it is positive. The finite energy content of the function is generally normalised to unity and is verified by: Z

∞

|ψ(y)|2 dy = 1

(2.12)

−∞

The final condition is known as the admissibility condition. This is required to ensure that the wavelet transform has an inverse and is expressed as, Z

∞ −∞

|Ψ(ω)|2 dy = Kψ < ∞ |ω|

where Ψ(ω) is the Fourier transform of the mother wavelet. 9

(2.13)

In this study the popular and purely real Mexican-hat wavelet is used (equation 2.14, figure 2.2). ψ(y) =

d2 − y 2 e 2 dy 2

(2.14)

Figure 2.2: The mexican hat mother wavelet. The wavelet transform of a 1-D function f (y) is defined as follows, 1 f˜(y, λ) = λ

µ

Z 0

f (y )ψ

∗

y0 − y λ

¶ dy 0

(2.15)

where ψ ∗ is the complex conjugate of the mother wavelet. The transform, f˜, is a function of y, a coordinate in 1-D space, and λ, the wavelet scale in the y-dimension. The wavelet transform is closely linked to the Fourier transform [5], which expresses a function in terms of trigonometric components, providing information on frequency content. The oscillations of the trigonometric components continue across all space, causing a delocalisation of the information contained in the original function. The advantage of the wavelet transform is that it maintains the locality of the information by using localised mother wavelets. Figure 2.3 shows an example one-dimensional function and its wavelet transform. The wavelet transform is equally applicable in higher dimensions using higher dimension mother wavelets.

10

(a)

(b)

Figure 2.3: a) An example function and b) plot of the wavelet transform coefficients as a function of space and scale. 2.3.2

The global wavelet energy spectrum

The energy associated with a given scale is described by the global wavelet energy spectrum, which is defined as follows. 1 G(λ) = λ

Z f˜(y, λ)2 dy

(2.16)

The spectrum expresses the distribution of a function’s energy across scale space, providing a description of the scale composition of the function. The wavelet transform contains a large quantity of information. The advantage of the global wavelet energy spectrum is that it condenses the information into a form that is often more appropriate for direct comparison.

11

2.4 2.4.1

Fractals Introduction

The term fractal refers to a geometric shape that cannot be described by traditional Euclidean concepts. The statistics of a standard geometry may be easily quantified; for example, the perimeter of a circle is a finite quantity and can be described as a function of its radius. Fractal geometries do not display such properties. The Koch flake, shown in figure 2.4, is an example of a fractal.

Figure 2.4: The Koch flake. The Koch flake is a mathematical fractal comprising infinite repeating scales. Every edge is interrupted by further repeating detail. The length of the perimeter of the Koch flake is thus infinite. The area, however, converges to a finite value. Fractals occur frequently in nature. Unlike mathematical fractals such as the Koch flake, naturally occurring fractals contain a finite number of scales, eventually becoming smooth at small scales. Turbulence exhibits notable fractal behaviour pertaining to the dissipation of energy through the scales of the energy cascade. The range of scales of the turbulent eddies is analogous to the repeating scales which characterise fractal geometry. In the study of turbulence, fractals are useful as standard Euclidean concepts are often inadequate [7].

12

(a)

(b)

(c)

Figure 2.5 Fractals in nature - a) a fern leaf, b) the coastline of Norway, c) the human lung.1 2.4.2

Fractal dimension

Dimensionality is a well known concept of Euclidean geometry and may be extrapolated to fractals. Fractal dimension provides a measure of fragmentation which is useful in the study of flow interfaces. Standard geometric shapes fall into four categories. Category

Dimension, D

point

0

line

1

area

2

volume

3

From these categories, it follows that an object of unit size and dimension, D, may be covered by N similar objects scaled by factor λ, according to relationship 2.17.

N = λ−D

(2.17)

For example, a 2-D square is covered by 4 similar squares scaled by a half, while a 3-D cube is filled by 8 similar cubes scaled by a half (see figure 2.6). 1

Images are re-used under the GNU Free Documentation License.

13

Figure 2.6: A demonstration of self-similar object coverage in 2- and 3-D. By rearrangement, the dimension of the object may be expressed as follows.

D=−

log N log λ

(2.18)

Application of this relationship to fractal objects demonstrates an interesting property. Cantor dust is one of the simplest fractal examples. Starting with a one-dimensional line of unit length, with each successive generation the middle third of every contiguous line is removed resulting in lines of one third the length as demonstrated in figure 2.7.

Figure 2.7: Cantor dust. For Cantor dust, N and λ can be be expressed as a function of the generation n as in 2.19. Generation, n N

λ

0

1

1

1

2

1 3

2

4

1 9

N = 2n , λ = 3−n

(2.19)

A priori, the fractal dimension of Cantor dust is less than one as subtractions are performed to a line of initially one dimension. Through application of relationship 2.18, the

14

fractal dimension of Cantor dust is found to be 0.63.

DCantor = −

log 2 n log 2 = ' 0.63 −n log 3 log 3

(2.20)

This dimension indicates the extent to which the Cantor dust fills one-dimensional space. The relationship between the dimension of a fractal object embedded in q-dimensional space to the dimension calculated from a p-dimensional cross-section can be shown to follow relationship 2.21 if the section is “independent of the fractal itself ”[7], holding for any section of an isotropic fractal.

Dq = Dp + (q − p)

(2.21)

For example, the true dimension of a 3-D isotropic flow field can be related to the dimension of a 2-D section, just as the dimension of a square can be related to the dimension of a cube; by the relationship D3 = D2 + 1. The fractal dimension of the mixing interface provides a measurement of folding - how much the interface fills space. A dimension close to 2 of an object embedded in 2-D space indicates that the object nearly fills that space. If a dimension of 2 is calculated for a 2-D slice of an object embedded in 3-D space, for an iso-tropic fractal, the object has a dimension of 3; indicating that it also fills 3-D space. Applied the species interface in a mixing field, it has relevance to measuring mixing.

15

Figure 2.8: A box-counting example. 2.4.3

Box-counting

Box-counting is a numerical method used to calculate the fractal dimension of an object by exploiting equation 2.18. The object of interest is covered with an array of boxes and those which contain a portion of the object are counted. With each successive iteration the size of the boxes is reduced. The number of occupied boxes, N , is plotted as a function of the reciprocal of the box edge length, λ1 . On a log-log scale, a constant gradient indicates a continuous distribution of scales and is equal to the fractal dimension.

D=

log N log N 1 = − log λ log λ

(2.22)

Figure 2.8 is an illustration of the box-counting method on the black/white interface of a binary map.

16

3

Procedure

A flow of marker fluid through fractal and non-fractal grids is observed using digital photography. The resulting flow patterns are analysed using the box-counting, wavelet and statistical methods discussed, to examine the effect of the grids on the flow in terms of fractal behaviour, scale development and mixing.

3.1 3.1.1

Experimental procedure Apparatus

An experimental apparatus has been constructed comprising a perspex channel, into which a grid can be fixed. A coloured marker fluid is introduced into a main-flow of water at the entrance to the channel. Immediately downstream of the grid, the flow is photographed in a viewing section using a tripod mounted digital camera (see figure 3.2). The flow channel has a rectangular internal cross-section measuring 162 mm × 20 mm. The upstream section is 1.25 m in length. Figure 3.1 shows a schematic of the apparatus. 3.1.2

Flow conditions

Flow through the apparatus is supplied by the reservoir and is driven by gravity. The flow rate must be limited to ensure an established flow before the grid. For a channel flow, the Reynolds number is expressed as,

Re =

ρudh µ

(3.1)

where dh is the hydraulic diameter and is defined as,

dh =

4A P

where A is the flow area and P is the perimeter of the flow area. 17

(3.2)

18 Figure 3.1: The experimental set-up.

For the purposes of calculation the flow is considered laminar, thus the relationship between dh and the development length, Le , required for flow to establish is as follows. Le ≈ 0.06Re dh

(3.3)

For the available length in the upstream channel, Le = 1.25 m, a maximum Reynolds number of Remax ≈ 585 is found. This is the around the upper limit for laminar flow and corresponds to a flow rate of 3.2 l min−1 (∼ = 1.6 cms−1 mean velocity in the channel). To meet the criterion of developed flow at the grid, a flow rate of 3 l min−1 is maintained. While the upstream flow is laminar, turbulence is generated by the grid. 3.1.3

Image capture

Compressed JPEG images are taken at the viewing section using a tripod mounted digital SLR2 camera as shown in figure 3.2. The digital camera specifications are included in appendix B.1. Images are captured for a period of around 200 s at the fastest rate obtainable with the camera. The capture rate slows from an initial 3 frames per second to one frame every 4 seconds due to time-to-disk limitations. A remote trigger is employed to avoid movement during image capture. All images are taken with a shutter time of 1/49 s, an ISO speed of 100 and an aperture value of 4.97 EV (f/5.6). Lighting is provided from below using a white laminated sheet as a diffuse background. The white balance is calibrated using an image taken with no marker species in the flow. Certain artefacts must be avoided to reduce interference with the image analysis. Following each experiment, the viewing section is cleaned to remove debris and dye deposits. As far as possible bubbles are removed from the viewing section prior to the experiment. A polarising filter and photographic umbrella are used to reduce reflections 2

Single Lens Reflex

19

from surrounding light sources. During analysis the images are assumed to represent a two-dimensional slice of the concentration field. Furthermore, the assumption is made that the image represents an orthographic view. One sequence of images is captured per grid (approximately 70 images). 3.1.4

Dye injection

Water soluble paint is used in the marker fluid. The paint is pre-mixed with water with a volume ratio of 1:150. The marker fluid is introduced into the main flow by direct injection using a syphon feed. The driving head of marker fluid, with respect to the level of the main reservoir, is fixed across all experiments. The resulting flow rate of dye is close to constant at 1.25 l min−1 . The entire dye tank capacity of 1.5 l is released during an experiment. Blue dye is used as it is easily distinguishable from the yellow tinted indoor reflections that are most prominent in the laboratory. This is to allow the option of employing a hue filter to the image to remove associated artefacts. A number of injection methods have been considered. Previous work by Brishni (2008) [15] was carried out using a submersible pump. This method presents a number of disadvantages: In order to submerge the pump, a large quantity of marker fluid is required and a significant quantity cannot be delivered. Furthermore, with the pump deactivated, a syphoning effect leads either to continued dye delivery or draw-back of the main flow into the dye tank. A syringe pump system offers accurate flow control and minimal wastage, but limits the volume of fluid that can be delivered. An expensive compromise is offered by a peristaltic pump. The syphon feed is inexpensive and delivers a near-steady flow, flexibility in the quantity of dye to be delivered, and minimal wastage.

20

Figure 3.2: The image capture set-up.

Figure 3.3: Control grid 1 in place with the viewing section detached.

21

3.1.5

The grids

A total of eight perspex grids are tested in four pairs; the fractal Cantor grid and its corresponding control grid, which has the same flow area but a non-fractal geometry. The four pairs correspond to extruded Cantor dust generations one to four (see section 2.4.2). Figure 3.4 shows control grids one to four and the fourth generation Cantor grid with the relevant dimensions. Figure 3.3 shows a grid in place, with the viewing section detached. Due to the constraints of the material properties and cutting tools, the minimum distance between flow areas is 2 mm. This limits the fractal grids to the fourth generation and constrains the shape of the control grids. The control grid geometries employed are arbitrary and are not unique solutions. Grid anatomy comprises flow areas, through which fluid may pass, and blockage areas which are impenetrable. The blockage areas lead to blockage regions in the immediate downstream flow.

3.2

Numerical analysis

The image sequences are prepared using Linux/Unix Bash scripts and Matlab and all image analysis is undertaken using Matlab. The University of Sheffield Iceberg server and Sun Grid Engine are employed where more intensive computing power is required. Prior to analysis, the relevant data is extracted from the digital images. Figure 3.6 shows the definition of the x and y axes used for the analysis of the digital images. In order to identify the coloured dye against the white background, the image encoding is converted from standard RGB (Red, Green, Blue) colour space to HSV (Hue, Saturation and Value) colour space (see [14, 13]). A map of saturation, S(m, n) ∈ [0, 1], is created, where m and n are the pixel coordinates. S(m, n) is a digital representation of the concentration field, C(x, y). The marker fluid shows high saturation values while the white background exhibits a near zero value. The saturation value is assumed to be

22

Figure 3.4: Dimensioned drawings of control grids 1, 2, 3 & 4, and the fourth generation Cantor grid. 23

indicative of dye concentration and the relationship is assumed linear. The conversion of digital images to saturation maps provides a clear representation of the concentration field, independent of shading effects. Figure 3.5 shows a comparison of the value map and saturation map methods and clearly shows the improved clarity and independence of lighting effects obtained using the saturation method. It is important to note that the saturation map may only store data in the range [0, 1]. To the fullest possible extent, the camera settings have been selected such that the full range of concentration values is captured within this range.

(a)

(b)

Figure 3.5: The concentration field represented by a) an inverted map of pixel value, and b) a map of pixel saturation. The saturation map, S(m, n), is analysed directly using the wavelet method. For box-counting, a binary map π(m, n) is created using a chosen threshold value, S 0 . The pre-processing procedure is illustrated in figure B.1.

π(m, n) =

   1 if S(m, n) ≥ S 0   0 if S(m, n) < S 0

(3.4)

The turbulent field is non-isotropic and three-dimensional; however, all analysis is carried out on a two-dimensional sample, which is assumed to be representative of the full field. To allow analysis of the development of the flow over time, the image capture time is 24

Figure 3.6: Axis definitions for the digital images extracted from the image meta-data, saved to file and read by the analysis algorithms. A Linux BASH script is employed to prepare the image sequences for Matlab. Using the ImageMagick library the images are rotated to compensate for any skew in the image capture angle and then cropped to a 15 × 20 cm region in the centre of the viewing area, using reference markings on the background. The image capture time is extracted from the meta-data using the EXIFTOOL library. Time t = 0 is defined as the point at which the instantaneous threshold, S 0 = 0.1. Separate Matlab codes are employed to run batch processing of all images in a sequence. Thresholding The choice of a threshold value is very important as it has a significant effect on the results of the box-counting. Two methods of thresholding are employed. For a given image, an optimum threshold, SO 0 , may be calculated. For this investigation, the Otsu method is used (Otsu 1979 [6]). A fixed concentration iso-level is also tracked. The two methods are as follows: Instantaneous threshold method Each of the images in a sequence is converted to a binary map πO (m, n) based on its Otsu threshold. Fixed threshold method The image sequence is first analysed to obtain the instantaneous Otsu threshold as a function of time. Each image is converted to a binary map, πm (m, n), based on the median threshold, SO 0 , of the entire sequence.

25

3.2.1

Box-counting

A box-counting algorithm has been developed to calculate the coordinates of the boundary between the black and white regions of the binary map π(m, n). An initial box edge length of

lB , 10

where lB is the largest edge of the image, is chosen

to ensure statistical validity. The covering box size, λ, decreases with each iteration, n, according to the relationship λ =

lB . n1.2

A minimum box size of 3 x 3 pixels is imposed to

avoid sampling beyond the image resolution. The algorithm examines each boundary point in turn to find the box into which it falls. The box is identified using a simple search algorithm. Many existing box-counting codes divide the image into 2n boxes. The result is that the boundary between two boxes always fall in the same place. The algorithm written for this investigation divides the image into boxes of size

lB . n1.2

The boundaries therefore

fall in different places with each iteration, avoiding any possible aliasing effects that may occur. 3.2.2

Wavelet analysis

A wavelet method has been developed, based on the global wavelet energy spectrum (see section 2.3). The spectrum of a one-dimensional section of the saturation map is calculated at given distance downstream from the grid. This section is taken in the cross-flow y-axis. The global wavelet energy spectrum of S(m0 , n) (where m0 is fixed) has a purely qualitative significance, and no direct physical meaning. 3.2.3

Statistical analysis

A separate code is employed to calculate the coefficient of variance and the exposure of each of the saturation maps in a sequence.

26

4

Results

Some example images are presented, showing the development of the flow through the first generation control grid. Subsequently the statistical, wavelet and box-counting analyses are presented for each equivalent pair of grids and are briefly described.

(a)

(b)

(c)

(d)

Figure 4.1: The saturation map of the flow through the first generation control grid at t = 2, 30, 70 and 120 s.

(a)

(b)

(c)

(d)

Figure 4.2: The binary map of the flow through the first generation control grid at t = 2, 30, 70 and 120 s, using the fixed value threshold method (SO 0 = 0.4941).

(a)

(b)

(c)

(d)

Figure 4.3: The binary map of the flow through the first generation control grid at t = 2, 30, 70 and 120 s, for the instantaneous Otsu threshold (SO0 = 0.0510, 0.5098, 0.8902, 0.7608).

27

4.1

Statistical properties

(a) Grid 1

(b) Grid 2

(c) Grid 3

(d) Grid 4

Figure 4.4: Exposure, E, Coefficient of Variance, CoV, and Otsu threshold as a function of time for equivalent grid pairs. At time t ' 0, the exposure is close to zero as the marker species is not present and the saturation map indicates a uniform zero concentration field (neglecting the effect of artefacts). Exposure rises steadily as marker enters the image, reaching a maximum as the dye begins to dominate the flow. The exposure drops as coverage increases, reducing the interface size and sharpness. As the marker injection ceases, the pattern repeats as the colourless flow returns. The coefficient of variance is initially high, as the species are well separated - dropping as mixing occurs and the intensity of segregation reduces.

28

4.2

Global wavelet energy spectrum

(a) Grid 1

(b) Grid 2

(c) Grid 3

(d) Grid 4

Figure 4.5: Global wavelet energy spectrum of a slice perpendicular to the flow, at x ' 10 cm from the downstream image edge for equivalent grid pairs. Figure 4.5 shows the global wavelet energy spectrum of 1-D cross-flow slices taken 10 cm from the downstream image edge. The slices are taken from a single image from each sequence. The time criterion for the single sample image is SO0 = SO0 . This occurs at approximately t = 20 s for all sequences. The peak energy generally occurs at a scale of 7 . λ . 8 cm. The first iteration peaks at a lower value, reflecting the entrainment of the flow into a single central flow. The fractal grids exhibit a step at lower scales.

29

4.3

Fractal dimension

(a) Grid 1

(b) Grid 2

(c) Grid 3

(d) Grid 4

Figure 4.6: Fractal dimension of the median threshold interface for equivalent grid pairs (see figure 4.2 for example images). Fractal dimension, D, calculated at the median threshold value, is close to zero at time t → 0. The fractal dimension rises sharply as the dye enters the image, reaching a peak. As the above-threshold region consumes the image, the fixed interface moves to the edges and the dimension tends towards one. The presence of artefacts in the image cause some interference, appearing as small above-threshold regions. This effect is negligible compared to the regions representing the dye and only manifests itself for a few seconds.

30

(a)

(b)

Figure 4.7: a) The correlation coefficient of the box-counting slope for control grid 1 between t = 0 and 20 s, b) comparison of fractal dimension of the flow through the first generation grids, calculated using the instantaneous threshold method. The correlation coefficient of the box-counting slope,

logN 1 , log λ

indicates its linearity, where

1 is a perfect linear gradient and 0 shows no correlation. Figure 4.7 a) is an example correlation taken from control grid 1. Linearity is high and from around 10 s maintains a value of 1. This is an indication of the distribution of scales and the applicability the box-counting method. Figure 4.7 b) shows the fractal dimension of the flow through the first generation control and fractal Cantor grid, calculated using the instantaneous optimum threshold as discussed in section 3.2. The CPU time required to complete this calculation for higher generations is significant (> 200 hours per grid) due to very large interfaces. As such only the first generation grids have been analysed using this method.

31

5

Discussion

Before the results are discussed and interpreted, it is useful to establish some definitions regarding the development of the concentration field over the sampling time window and within the sampling frame (the image). These definitions are constructed from general observation and related to the analytical results. The significance of the measurement data is considered and a comparison of the grids is made. Finally, the assumptions and limitations of the method are reviewed. Where these limitations are relevant to the understanding of the results, they will be mentioned as appropriate.

5.1

Four stages of development

For the purposes of description, the term flow describes the velocity vector field and concentration field is a scalar field describing species distribution. The marker fluid is inert and, being composed of water soluble paint in water, exhibits properties close to that of the main flow; therefore, the problem is classed as level 1 mixing - the fluid mechanics are independent of the mixing process. Although the mixing problem is unsteady and anisotropic, the main flow is stationary. Marker fluid is introduced into the downstream main flow and quickly mixes in a turbulent regime, becoming laminar before reaching the grid.

Figure 5.1: Image of the flow regime for the first generation fractal grid (manipulated for clarity). 1) Homogeneous incident concentration, 2) grid flow area, 3) grid blockage area, 4) blockage region, 5) flow region. Figure 5.1 shows the flow regime around a grid. The incident flow (immediately 32

upstream of the grid) is by calculation and observation, laminar. The fully established incident concentration field is close to homogeneous. This homogeneous state is interrupted by the grid. An unsteady mixing regime is established whereby the entrained flow develops towards re-establishing a homogeneous concentration. The entrained flow (immediately downstream of the grid) is highly dependant on the grid geometry and conditions range from weakly transitional to turbulent. The number of discrete flow areas in a grid determines the number of individual flows that comprise the entire entrained flow. These discretised flows shall be labelled subsidiary flows. The development of the concentration field may be described in four chronological stages, during which different mechanisms may be observed. Entry Marker fluid begins to enter the frame, entrained by the grid. During this stage the entrainment mechanism is easily observed. Entry is characterised by a peak in the intensity of segregation (as described by the coefficient of variance) and a steady increase in the exposure. During entry a sharp rise in the fractal dimension of the fixed value interface occurs, from zero, indicating no interface (neglecting artefacts in the saturation map) to a value between 1 and 2. Mixing Subsequently dye bridges the entire frame. The kinematic stirring mechanism can be observed as the intermediate turbulent scales distribute the dye. The dye interface, however, remains well defined as the visible effects of molecular mixing are thus far negligible. The result of a well defined interface is a peak in exposure. The intensity of segregation begins to fall as the marker region grows; the mean saturation value increases and thus the variance of the saturation map reduces. As the mixing stage begins, the fractal dimension of the fixed value interface peaks. Uniformity As the molecular mixing mechanisms act, the mixing interface becomes less defined and the concentration field approaches a homogeneous state. In the cases where a fully homogeneous state is achieved, it can be considered that full mixing has 33

occurred. Uniformity is characterised by a plateau in the exposure and intensity of segregation. In a perfectly uniform saturation map, the exposure and coefficient of variance are equal to zero. As the fixed value interface migrates to the edges of the frame, its fractal dimension, D → 1. Clear-down A finite quantity of dye is released and eventually the incident concentration diminishes and regions of low concentration enter the frame. This stage can be considered a similar repetition of the earlier stages.

5.2

Statistical properties

The significance of the development of exposure and intensity of segregation is discussed. The scales of segregation are discussed separately in section 5.3. As expected, the intensity of segregation provides a good indication of mixing with reference to uniformity. Where a homogeneous concentration field is reached, the coefficient of variance reaches a value very close to zero. It is important to mention the discrepancy between the coefficient of variance and the intensity of segregation that it describes. At time t = 0, there exists a homogeneous concentration field containing no marker species, thus IoSt=0 = 0; however, in practice, the saturation map is sensitive to artefacts present on the channel walls, such as scratches, bubbles or debris, which may appear as distinct regions of a second species and CoVt=0 > 0. Exposure responds to two specific factors which serve to increase the length of the dye interface. The discretisation of the flow by entrainment at the grid leads directly to a larger interface. For example, the first generation Cantor grid entrains flow into a single subsidiary flow with two interface origins (one at either extremity of the flow area). The equivalent control grid entrains three subsidiary flows with a total of six interface origins. The kinematic stirring mechanism, induced by recirculation in the blocked regions, also increases the complexity of the interface, increasing its length. 34

Molecular and viscous diffusion serve to reduce the sharpness of the interface. This manifests a reduction in exposure. High exposure coincides with clarity of the kinematic stirring mechanism; however, it is an indication of the concentration gradient and according to Fick’s law must therefore be linked to the mass flux attributed to the diffusion processes.

5.3

Wavelet analysis - the scales of segregation

The one-dimensional global wavelet energy spectra are intended to provide an indication of the scale composition of the saturation map and thus an insight into the scales of segregation of the concentration field. In order to reduce the information sufficiently for comparison, the spectra are fixed in time and space. Difficulty arises in specifying an equivalent time location for each of the spectra. As the problem is unsteady and involves a large number of variables, it is difficult to define an equivalent time for each grid. An attempt is made to select equivalent time locations according to the saturation threshold; however, this method does not necessarily produce comparable results. The results do not show any particular patterns that suggest comparability between data sets. The wavelet method may be more appropriately applied to a stationary mixing problem. Nevertheless, the results allude to some interesting phenomena which may be confirmed by a stationary experiment. This is discussed as part of the grid comparison (section 5.5). The possibility of estimating the length of the Kolmogorov scale is shown in section 2.1.1; however, it is necessary to know the size of the integral scales. It may be useful to compare these calculated scales with the global wavelet energy spectra. In general, flow velocity measurements are required to measure the integral scales, although they can be observed in some, but not all, of the images. This is an unfortunate limitation of the non-invasive techniques used and the acquisition of only scalar data.

35

5.4

Fractal dimension

The two thresholding methods discussed in section 3.2 lead to very different box-counting results. Instantaneous threshold method The instantaneous threshold interface is very large for a near homogeneous saturation map. The CPU time required to compute the development of the fractal dimension of such a large interface exceeds 200 hours per grid. The first generation grid flows do not reach a homogeneous state and can be analysed using the instantaneous threshold method, with a CPU time of around 30 hours. The absence of results for higher generation grids using this method prevents comparison and interpretation. Fixed threshold method The development of the fixed interface dimension illustrates well the four stages of mixing and displays similar attributes to the exposure value. The fixed value interface method imposes a binary definition of mixing; all sub-threshold regions are defined as unmixed while all above-threshold regions are mixed. This simplistic view limits the significance of fractal dimension as a measurement of mixing. One significant result is that the fractal dimension of the fixed value interface peaks consistently at 1.35 . D . 1.4. During their investigations into the fractal nature of turbulence, Sreenivasan and Meneveau (1986) consistently found the dimension of the turbulent/non-turbulent interface to be in the region of 1.3 to 1.4 for a two-dimensional slice of a fully established isotropic turbulent flow field. It is also worth noting that this peak closely coincides with the median sequence threshold value - the peak occurs when the fixed threshold is equal to the instantaneous optimum threshold. This observation implies that the optimum threshold value indicates well the turbulent/non-turbulent interface. It is possible that some information may be acquired from the slope of the fixed value interface dimension as it rises from an initial value of zero. However, low time resolution

36

makes comparison difficult. Furthermore, although no detailed study has been conducted, the linearity of the box-counting calculations appears to be relatively low during the entry stage, implying that a fractal dimension may not be applicable. A more in-depth study of the correlation coefficients of the box-counting results would be required to assess the applicability of a fractal dimension at this stage.

5.5

Grid comparison

The comparison of the statistical properties of the grids shows that higher generation grids behave similarly to their equivalent control counterparts. Furthermore, the results show that faster mixing occurs in higher generation grids; however, this result is biased by the nature of the problem. The incident concentration field is homogeneous and since the higher generation grids have a minimal effect on the flow, the entrained flow returns very quickly to this state. The first generation grids have large blockage areas at the extremities. The effect of this geometry is large blockage regions in the entrained flow which inhibits the development of a fully homogeneous downstream concentration. This exacerbates the bias towards apparently faster mixing in higher generation grids. With reference to the various measurements taken, some marked effects can be observed, especially in the first generation grids. The exposure of the first generation control grid reaches a significantly higher value (E ' 0.024) than the other grids. This implies improved kinematic stirring and therefore improved molecular diffusion. The grid is defined by flow discretisation and large blockage areas which induce strong turbulent recirculation. The first generation fractal grid, containing blockage regions but no flow discretisation, reaches a significantly lower peak exposure of around E = 0.017. Higher generation grids with high discretisation but very small turbulent induction due to small blockage areas, peak at up to E = 0.02. Observation of the saturation maps indicate that the flow around the blockage areas of the third and fourth generation grids is close to laminar. This observation suggests that 37

there may be balance to achieve between discretisation and blockage area size. Further investigation would be required to test this hypothesis. Further observation of the flow field suggests that the most important factor in the induction of turbulence is the geometry of the blockage areas rather than flow areas. The Cantor grid model used in this investigation is such that, while the flow areas are of different sizes the blockage areas are of the same size. A more significant effect may be imposed on the flow by choosing an inverted model whereby the blockage areas are varied and the flow areas constant. This may effect a greater impact on the flow development. An inverted Cantor model would lead to very small flow areas. Different fractals could be explored to increase the flow area. Although comparison of the global wavelet energy spectra is difficult, there are some apparent differences in the scale distributions in the saturation maps of the fractal and control grids. The spectra show a small step at the small scales which appears to be exaggerated in the fractal results. A stationary experiment would allow this to be investigated more thouroughly by eliminating variables for better comparison. Interestingly, while the flow through the fourth generation grids are, by observation, laminar, the fractal dimension nevertheless peaks at D ' 1.35, as would be expected for a fully turbulent interface. This tends to contradict the idea that turbulence is a special case in terms of fractal analysis of the interface.

5.6

Method review

It is important to reflect upon the assumptions made during the investigation and the possible effects on the results. These assumptions are discussed, followed by a review of the limitations of the experimental procedure. There is a disparity between the method of image capture and the subsequent image analysis which arises from the assumption that the saturation maps generated from the digital images represent an orthogonal view of a two-dimensional slice of the concentration field. In reality they represent a perspective projection - they can be considered as many 38

slices superimposed upon one-another - the detail in a slice of fluid is masked by the layers around it. Box-counting calculations are affected by the masking of internal boundaries and possible smoothing effects on the external boundaries. Exposure results are affected by the cumulative effect of the projection which may cause the flow to appear more homogeneous - exposure may be under-calculated. A closer approximation to a twodimensional slice could be achieved using a laser light sheet, which would also reduce the effect of artefacts present on the channel walls. Bubbles in the viewing section are also a common source of artefacts as they refract light. Where possible most bubbles are removed from the viewing section but this is not possible in all cases. Bubbles close to the grid on the upstream side may also affect the flow conditions. A vertically orientated flow channel could alleviate such problems. The application of single statistical scalars (dimension, coefficient of variance and exposure) to the entire frame is incompatible with the nature of the problem, which is anisotropic and inhomogeneous. It may be more appropriate to divide the frame into zones, or take sample slices normal to the mean flow. Application of two-dimensional wavelet analysis with a zoned system may also be a useful method of analysis. The application of a 2-D sampling regime is also in conflict with the 3-D nature of turbulence. Certain difficulties arise in the definition of a time origin. An approximation is made with relation to the threshold value; time t = 0 is defined as the point at which the instantaneous threshold, SO 0 = 0.1. The sampling time resolution is also limited by the image capture equipment and inhibits the study of transient effects. Furthermore, the large sampling time required leads to a limitation of time resolution due to the practicalities of analysing large amounts of data. The time variable would be negated altogether by stationary flow conditions. Marker injection at the grid would create a stationary problem and allow better comparison of the grid effects. The relationship between the saturation map and the concentration field is assumed to be linear, although this has not been proven. Calibration could be performed on known homogeneous concentrations under the same photographic conditions in order to quantify 39

this link. A major limitation of the experiment is the restricted sampling frame. As the energy cascade develops in the direction of flow, small scale turbulence and molecular mixing effects occur downstream of the viewing section and are not present in the results. Analysis should be carried out on these downstream effects. It has also been observed that the flow is not consistently turbulent. Higher order grids display largely laminar properties. This is restricted by the development length imposed by the length of the upstream channel. An increased flow rate (with a longer upstream channel) or altered grid patterns could improve turbulent induction. A single set of control grids is used and the designs are not unique. Further designs should be tested.

40

6

Conclusion and Recommendations

An investigation has been carried out into turbulent mixing of a marker fluid in a flow through fractal Cantor grids using image analysis. Following previous work, the experimental procedure has been improved, particularly with regard to image capture techniques and the inclusion of control experiments using non-fractal grids. A set of analytical tools has been developed, improving existing box-counting methods and developing wavelet and statistical methods. Fractal cantor grids have been compared to non-fractal control grids of equivalent flow area. The problem consists of an incident flow of homogeneous concentration, entrained by a grid into a transient mixing field. The flow conditions in the mixing field range from weakly transitional/laminar to turbulent. Image sequences of the entrained flow have been converted to saturation maps which are assumed to represent an orthoganol 2-D slice of the marker species concentration field. The saturation map sequences have been analysed using image processing techniques. The intensity of segregation (represented by the coefficient of variance) and exposure have been calculated. An attempt has been made to examine the scales of segregation using global wavelet energy spectra of 1-D slices of the saturation map at a fixed time location. The fractal dimension of a fixed value saturation threshold has been calculated using the box-counting method. The following conclusions have emerged from the analyses and observation of the concentration field. ˆ The results suggest that exposure shows a strong positive response to turbulent

kinematic stirring and discretisation of the entrained flow. Both serve to increase the length of the interface, contributing to stretching, folding and break-up of the marker species. ˆ Exposure is linked to the concentration gradient and by Fick’s law is shown to be

a measure of the mass flux associated with molecular mixing. 41

ˆ Wavelet analysis alludes to increased energy at the small scales in the flow through

the fractal grids. ˆ The iso-value interface dimension peaks at 1.35 . D . 1.4, which is in agreement

with the findings of Sreenivasan and Meneveau (1986). This result is equally observed in the cases where the flow regime is observed to be near-laminar. Any special fractal properties for turbulent cases is not observed. ˆ Mixing may benefit from a balance between flow discretisation and kinematic stir-

ring induced by larger blockage areas at the grid. Further conclusions may be drawn regarding the improvements made to the experimental procedure. ˆ The use of a saturation map in lieu of a value map reduces the effects of non-

uniform light and shadow and more effectively identifies the marker fluid against the background. ˆ A polarised filter coupled with a photographic umbrella serves to render the effect

of reflections negligible. The results obtained from this investigation successfully highlight necessary improvements to the method. ˆ The application of single scalar measurements to the sampling frame is inappropri-

ate. A zoning scheme should be considered. ˆ A more comprehensive study of the applicability of the box-counting algorithm

should be conducted, examining the linearity of the box-counting results. ˆ A stationary problem could be created by injecting the marker species at the grid,

eliminating the time variable and allowing better comparison. Time-averaging could be employed. 42

ˆ A quantified relationship between the saturation map and concentration field should

be obtained. ˆ The sampling frame should be extended to capture the downstream flow. ˆ A longer upstream section should be used to allow higher flow rates. ˆ A wider set of fractal grids should be investigated including different fractal models,

inverted fractal models and varied designs of control grid. ˆ The effect on mixing of the size of the blockage areas and discretisation of the flow

areas should be investigated further. ˆ A laser light sheet, or similar method, should be employed to achieve a closer

approximation to a 2-D slice. ˆ Implementation of a vertical flow channel should be considered to mitigate bubble

accumulation.

43

References [A] Papers [1] Paul E. Dimotakis. “Some issues on turbulent mixing and turbulence”. In: GALCIT Report FM93-1a (Aug. 1998). [2] Paul E. Dimotakis. “Turbulent Mixing”. In: Annual Review Fluid Mechanics 37 (2005), pp. 329–356. [3] Marie Farge. “Wavelet transforms and their applications to turbulence”. In: Annual Review Fluid Mechanics (1992), pp. 395–457. [4] Alena Kukukova. “A new definition of mixing and segregation: Three dimensions of a key process variable”. In: Chem Eng Res Des (2009). doi: 10.1016/j.cherd. 2009.01.001. [5] F. Nicolleau and J. C. Vassilicos. “Wavelet analysis of wave motion”. In: Personal communication, F. Nicolleau, University of Sheffield (2000). [6] Nobuyuki Otsu. “A threshold selection method from gray-level histograms”. In: IEEE Transactions on systems, man, and cybernetics 9.1 (Jan. 1979), pp. 62–66. [7] K. R. Sreenivasan and C. Meneveau. “The fractal facets of turbulence”. In: Journal of Fluid Mechanics 173 (1986), pp. 357–386.

[B] Books [8] C. W. Fetter. Contaminant Hydrogeology. 2nd. Prentice Hall, 1999. Chap. 3. [9] Jean Mathieu and Julian Scott. An Introduction to Turbulence. Cambridge University Press, 2000. [10] J. M. Ottino. The Kinematics of Mixing: Stretching, Chaos, and Transport. Cambridge University Press, 1987.

44

[11] E. L. Paul, V. Atiemo-Obeng and S. M. Kresta. Handbook of Industrial Mixing: Science and Practice. John Wiley & Sons Inc., 2003. [12] Stephen B. Pope. Turbulent Flow. Cambridge University Press, 2000. [13] Milan Sonka, Vaclav Hlavac and Roger Boyle. Image Processing, Analysis, and Machine Vision. Third. Thomson, 2008.

[C] Online [14] Richard E Riley. A comparison of methods for identifying a dyed fluid in an image. www.reriley.co.uk/idfluid.php. 2008. url: www.reriley.co.uk/idfluid.php.

[D] Masters thesis [15] Brishni Mukhopadhyay. “Image analysis of flow through a Cantor grid”. MA thesis. The University of Sheffield, 2008.

45

Appendices A

Image sequences

The following pages contain the saturation map sequences for each grid.

46

Figure A.1: Saturation map sequence for control grid 1.

47

Figure A.2: Saturation map sequence for fractal grid 1.

48

Figure A.3: Saturation map sequence for control grid 2.

49

Figure A.4: Saturation map sequence for fractal grid 2.

50

Figure A.5: Saturation map sequence for control grid 3.

51

Figure A.6: Saturation map sequence for fractal grid 3.

52

Figure A.7: Saturation map sequence for control grid 4.

53

Figure A.8: Saturation map sequence for fractal grid 4.

54

B

Ancilliary information Table B.1: Digital camera specifications Manufacturer: Model: Resolution: Lens:

Canon EOS 300D 6.2 Mega pixels Canon Zoom EF-S 58 mm Ø 18-55 mm 1:3.5 - 5.6

55

Figure B.1: The image processing data flow.

56

C

Algorithms

This appendix contains the algorithms used for image processing. Slightly modified versions of these codes are used on the University of Sheffield Iceberg Server where augmented computing power is required. All code is the work of the author.

C.1

Time meta-data extraction

The following Linux BASH code reads the meta data from each image in a sequence and writes to file. # Dialogue to choose file srcdir=$(zenity --file-selection \ --directory --title "Choose containing directory...")

# Make a directory for processed data files if one does not already exist mkdir -p "$srcdir"/prep/

# Extract creation date and time and read into a temporary file exiftool -p ’$filename $CreateDate’ -q -f \ "$srcdir" >> "$srcdir"/prep/timedata.tmp

# Sort the data and write to file, delete tmp file sort -d "$srcdir"/prep/timedata.tmp >> "$srcdir"/prep/timedata.dat rm "$srcdir"/prep/timedata.tmp

# Rotate and crop... # Take parameters rotangle=$(zenity --entry --text "Enter CW rotation angle in degrees") 57

cropstr=$(zenity --entry --text "Enter crop string in the form \ wxh+x+y (i.e. size then offset) e.g. 1208x1312+403+198")

rsiz=$(zenity --entry --text "Enter output size (%)")

# Run mogrify commands an feed output to a progress dialogue box ( echo "# Copying images..." ; sleep 1 cp "$srcdir"/*.JPG "$srcdir"/prep/ echo "25" ; sleep 1 if [ "$rotangle" -ne "0" ] then echo "# Rotating..." ; sleep 1

mogrify -distort SRT "$rotangle" \ "$srcdir"/prep/*.JPG

echo "50" ; sleep 1 fi echo "# Cropping..." ; sleep 1 mogrify -crop "$cropstr" "$srcdir"/prep/*.JPG echo "75" ; sleep 1 echo "# Resizing..." ; sleep 1 mogrify -resize "$rsiz"% "$srcdir"/prep/*.JPG echo "100" ; sleep 1 echo "# Preparation complete. Click ’OK’." ; sleep 1 ) | 58

zenity --progress \ --title="Image preparation" \ --text="Initialising..." \ --percentage=0

59

C.2

Image conversion

The following Matlab code converts the digital image to a saturation map or binary map.

1

% Image to concentratuion map.

2

% Author: Richard E Riley 20/11/2008

3 4

% Converts a colour image to an intesity scale using saturatuion value to

5

% determine the concetration of a coloured dye.

6 7

% Call as in the following example:

8

% im2conc(image,conversion type,level)

where:

9 10

% is the image to analyse

11 12

% is the type of conversion to perform. Possible options are:

13

% 'bin' = binary conversion. This option should be followed by ; the

14

% greyscale threshold. a value of 0 uses the Otsu Method

15

% 'scale' = greyscale conversion where pixel intensity indicates

16

% dye concentration

17 18

function [image,olev] = im2conc(rgb image,conversion,level)

19 20

% Convert to HSV, to greyscale based on saturation value.

21

hsv image = rgb2hsv(rgb image);

22

s = hsv image(:,:,2);

23

if strcmp(conversion,'bin') % Compare instruction string

24

% Convert to binary using Otsu's method to find the threshold. if

25

% no threshold is specified.

26

if level == 0

27

level=graythresh(s);

28

disp(['Otsu Level = ' num2str(level)]);

29

end

60

30

image = im2bw(s,level);

31

olev = level;

32

elseif strcmp(conversion,'scale')

33

% Simply report greyscale image

34

image=s;

35

else disp('Please specify conversion type as bin or scale');

36 37

end

61

C.3

Box-counting

The following Matlab code peforms the box-counting calculation on a single image.

1

% A box−counting algorithm Author: Richard E Riley 26/02/2009

2 3

% Call as in the following example:

4

% boxcount(image,xdim,ydim,iterations,level)

5

% to analyse <xdim> & are the dimensions of the image in cm

6

% is the number of iterations to perform, this is limited to a

7

% 3x3 pixel box. is the saturation threshold (between 0 and 1). a

8

% value of 0 uses the Otsu's Method

where: is the image

9 10

% The box counting size decreases according to the number of iterations N

11

% with the following relationship: 1/ Nˆ−1.2

12 13

% Function returns [dimension,correlation] where dimension is the final

14

% dimension and correlation the correlation coefficient of the box−count

15

% result.

16 17

function [dimension,correlation] = boxcount mk2(rgb image,x real,y real, Nmax,level)

18 19

% Convert to HSV and to greyscale based on saturation value, and then to

20

% binary.

21

[bw,olev] = im2conc(rgb image,'bin',level);

22 23

% Combine and plot all boundaries

24

B = bwboundaries(bw);

% Detect boundaries

25 26 27

% Plot image and boundaries

62

28

figure(1); clf; imshow(rgb image); hold on;

29

for k=1:length(B)

30

boundary = B{k};

31

plot(boundary(:,2), boundary(:,1), 'g','LineWidth',2)

32

end

33 34

% Convert from cell data into a matrix of x−y coordinates describing

35

% all distinct boundaries.

36

B=cell2mat(B);

37

if B

38 39

% Scale coordinates to physical dimensions.

40

[b,a]=size(bw);

41

sfx=double(x real/a);

42

sfy=double(y real/b);

43

B real = [];

44

B real(:,2)=B(:,2).*sfx;

45

B real(:,1)=y real−(B(:,1).*sfy);

% Image origin is the top−left

whereas plot origin is bottom−left 46 47

% Here the bounding box is define as the size of the image. In fact, this is a

48

% legacy to older code.

49

xmin=0;

50

xmax=a*sfx;

51

ymin=0;

52

ymax=b*sfy;

53 54

% Plot boundary points on a real scale and display the bounding box and

55

% gives some information

56

figure(2); clf; plot(B real(:,2), B real(:,1),'b.','markersize',3);

57

axis([0 x real 0 y real]);

58

gbox(1) = rectangle('Position',[(xmin),(ymin),(xmax−xmin),(ymax−ymin)], 'EdgeColor','r');

63

59

hold on;

60 61

% Box−counting...

62 63

% Count boundary points

64

np=size(B);

65

np=np(1);

66

disp([num2str(np) 'boundary points']);

67 68

% Redefine origin to that of the bounding box

69

% Padded with the equivalent length of one pixel (sf−).

70

B real(:,2) = B real(:,2) − xmin + sfx;

71

B real(:,1) = B real(:,1) − ymin + sfy;

72 73

% Define the inital variables for box counting

74

BBx=xmax − xmin + sfx;

% Length of bounding box in x

75

BBy=ymax − ymin + sfy;

% Length of bounding box in y

76

BB=max([BBx,BBy]);

% Choose largest edge (thus all boxes will be

square) 77

fullbox=[];

78

bcount=[];

% Matrix to hold the grid references of all full boxes % Matrix to hold box size vs. the number of full boxes

79 80

% Start at box size smaller than bounding box so that enough boxes

81

% exist to give a good statistical average.

82

lbs=BB/10;

83

sf=1.2; % Box shrink factor

84 85

% Limit max iterations to give a minimum box size of 3x3 pixels.

86

if (lbs/Nmaxˆsf)<((max([sfx,sfy]))*3)

87

Nmax=floor((lbs/(3*max([sfx,sfy])))ˆ(1/sf));

88

disp(['Iterations limited to ' num2str(Nmax)]);

89

end

90 91

% Calculate number of calcs needed and create a progress dialogue. This

64

92

% provides only an estimation of progess. The main purpose is to

93

% reassure the user that the program is advancing.

94

calcs=np*Nmax;

95

prog = waitbar(0,'Please wait...');

96

itprog = waitbar(0,'Iteration progress...');

97 98

% For each box size step

99

for N=1:Nmax

100

% Define length of each box edge. Power increases reduction steps

101

% with each iteration.

102

lb=lbs/Nˆsf;

103 104

% For each boundary point, find which box it falls into.

105

for k=1:np

106 107

bx=1;

108

while B real(k,2) > (bx*lb) bx=bx+1;

109 110

end

111

by=1;

112

while B real(k,1) > (by*lb) by=by+1;

113

end

114 115 116

fullbox(k,:)=[bx,by]; % Read grid reference into array

117

waitbar((k/np),itprog,['Point ' num2str(k) ' of ' num2str( np)]);

118

end

119

fullbox=unique(fullbox, 'rows'); % Allow each full cell to appear only once

120

[f,g]=size(fullbox); % f is now equal to the number of full boxes

121

bcount(N,:)=[(1/lb),f];

65

waitbar((N*np)/calcs,prog,['Iteration ' num2str(N) ' of '

122

num2str(Nmax)]); % Update progress dialogue 123 124

% Plot the boxes which contain the boundary giving some

125

% reassurance of what the script is actually doing.

126

delete(gbox(:)); % Delete previous boxes

127

gbox=[];

128

for k=1:f

129

gbx=((fullbox(k,1)−1)*lb) + xmin;

130

gby=((fullbox(k,2)−1)*lb) + ymin;

131

gbox(k) = rectangle('Position',[gbx,gby,lb,lb],'EdgeColor', 'r'); end

132 133 134

% Use a single order poly fit and take first order coefficient

135

% to obtain fractal dimension

136

if N > 1

137

p = polyfit(log(bcount(1:N,1)),log(bcount(1:N,2)),1);

138

dimension = p(1);

139

dimlog(N,:)=[N,p(1)]; else

140 141

p(1)=0;

142

dimension = p(1); end

143 144 145

end

146

close(itprog);

147

close(prog);

148 149

% Final data processing

150 151

% Plot log−log of dimension as a function of box size

152

figure(3); clf; loglog(bcount(:,1),bcount(:,2));

66

grid on; xlabel('1/\epsilon (cmˆ−ˆ1)','fontsize',12); ylabel('Fractal

153

dimension','fontsize',12) 154 155

% Return final dimension (for ease of use as opposed to accessing the

156

% dimlog matrix.

157

dimension = p(1)

158 159

% Calculate the correlation coefficient of the box−count to indicate

160

% linearity and appicability of method to given image.

161

correlation=corrcoef(log(bcount))

162

corsiz = size(correlation)

163

if

corsiz(1) == 1 correlation=0

164 165

else

166

correlation=correlation(1,2)

167

end

168 169 170

% What to do if no boundary points exist.

171

else

172

dimension=0;

173

correlation=NaN; end

174 175 176

end

67

C.4

Box-counting batch code

The following Matlab code invokes the box-counting script for a sequence of images and includes time data.

1

function dimensions = bcbatch(xdim,ydim,Nmax,level)

2

% BCBATCH Executes box−counting on multiple images and appends time data.

3

% See boxcount.m for variable meanings.

4 5

% Show dialogue to select the images in the sequence.

6

[fname,pname] = uigetfile('*.jpg;*.JPG;*.bmp','Select containing folder','/ home/rriley/Documents/Work/Year 4/Final Project/MatLab/Data/',' multiselect','on')

7

fname = sort(fname);

8 9

% Read the time data and sort into and array. Number of data lines should

10

% equal the number of images in the sequence.

11

[files, dates, times] = textread([pname 'timedata.dat'], '%s %s %s');

12

timedata=horzcat(files, times);

13

timedata = sortrows(timedata,1);

14 15 16

dimensions=[];% Create array to hold the time and calculated data.

17 18

% For each image invoke boxcounting algorithm

19

for k=1:length(fname);

20

f=cell2mat(fname(k))

21

p=([pname f]);

22

rgb image=imread([pname f]);

23

[dimensions(k,2),correlations(k)] = boxcount mk2(rgb image,xdim,ydim, Nmax,level);

24

time=char(timedata(k,2));

25

dimensions(k,1)=datenum(time,'HH:MM:SS');

68

26

end

27 28

% Provide plots and save data to file.

29

dimensions(:,1)=(dimensions(:,1)−min(dimensions(:,1)))*24*60*60;

30

plot(dimensions(:,1),dimensions(:,2),'o');

31

xlabel('Time (s)','fontsize',12); ylabel('Fractal dimension','fontsize',12) ;

32

fid = fopen([pname 'dimdev.dat'], 'w');

33

data=dimensions';

34

fprintf(fid, '%1.4f\t%1.4f\n', data);

35

fclose(fid);

36 37

plot(dimensions(:,1),correlations(:),'o');

38

xlabel('Time (s)','fontsize',12); ylabel('Correlation Coefficient',' fontsize',12);

39

fid = fopen([pname 'correlation.dat'], 'w');

40

dimensions(:,1)

41

correlations

42

data=[dimensions(:,1),correlations'];

43

fprintf(fid, '%1.4f\t%1.4f\n', data);

44

fclose(fid);

69

C.5

Statistical analysis

The following Matlab code calculates the coefficient of variance and exposure of all images in a sequence and includes time data.

1

function [CoV,E] = mixdim(name)

2

% [CoV,E] = mixdim(name) calculates the Coefficient of Variance (CoV) and

3

% Exposure (E) of an image sequence. Data is output to tab delimited file

4

% as a function of time. Time stamps are read from the timedata.dat file

5

% included in the folder containing the image sequence. (timedata.dat is

6

% generated using a bash script). Output is saved to the same directory as

7

% the image files. Variable 'name' is the label to be attributed to the

8

% output files.

9 10

% Show dialogue to select the images in the sequence.

11

[fname,pname] = uigetfile('*.jpg;*.JPG;*.bmp','Select containing folder','/ home/rriley/Documents/Work/Year 4/Final Project/MatLab/Data/',' multiselect','on')

12

fname = sort(fname);

13 14

% Read the time data and sort into and array. Number of data lines should

15

% equal the number of images in the sequence.

16

[files, dates, times] = textread([pname 'timedata.dat'], '%s %s %s');

17

timedata=horzcat(files, times);

18

timedata = sortrows(timedata,1);

19 20

dimensions=[]; % Create array to hold the time and calculated data.

21 22

% For each image perform calculations.

23

for k=1:length(fname);

24

f=cell2mat(fname(k))

25

p=([pname f]);

26

rgb image=imread([pname f]);

70

27

conc = im2conc(rgb image,'scale',0); % Convert image to saturation map

28 29

% Intensity of Segregation calculation

30

mconc=mean(mean(conc));

31

I = ((mconc − conc).ˆ2)./(mconc*(1−mconc));

32

dimensions(k,2) = sqrt(mean(mean(((conc−mconc)./mconc).ˆ2)));

33 34

% Exposure calculation

35

[a,b] = size(conc);

36 37

% Begin by adding one pixel padding to the saturation map. Padding

38

% is a copy of the border pixels to achieve zero difference between

39

% image edges.

40

conc2 = conc;

41

conc2=cat(1,conc2(1,:),conc2);

42

conc2=cat(1,conc2(a+1,:),conc2);

43

conc2=circshift(conc2,[−1,0]);

44 45

conc2=cat(2,conc2(:,1),conc2);

46

conc2=cat(2,conc2(:,b+1),conc2);

47

conc2=circshift(conc2,[0,−1]);

48

[a,b] = size(conc2);

49 50

% Reset value and begin exposure calculation

51

d=0;

52

for m=2:a−1

53

for n=2:b−1

54

p=[];

55

p(1) = abs(conc2(m,n) − conc2(m−1,n));

56

p(2) = abs(conc2(m,n) − conc2(m+1,n));

57

p(3) = abs(conc2(m,n) − conc2(m,n+1));

58

p(4) = abs(conc2(m,n) − conc2(m,n−1));

59

d=d+0.5*sum(p);

60

end

71

61

end

62

dimensions(k,3)=d/((a−2)*(b−2)); % Normalise and store in matrix.

63 64

% Send time to matrix.

65

time=char(timedata(k,2));

66

dimensions(k,1)=datenum(time,'HH:MM:SS');

67

end

68 69

% Provide plots and save data to file.

70

dimensions(:,1)=(dimensions(:,1)−min(dimensions(:,1)))*24*60*60;

71

subplot(2,1,1); plot(dimensions(:,1),dimensions(:,2),'o');

72

xlabel('Time (s)','fontsize',12); ylabel('CoV','fontsize',12); hold on;

73

subplot(2,1,2); plot(dimensions(:,1),dimensions(:,3),'x');

74

xlabel('Time (s)','fontsize',12); ylabel('Exposure','fontsize',12);

75

fid = fopen([pname name '−mixdim.dat'], 'w');

76

data=dimensions'

77

fprintf(fid, '%1.4f\t%1.4f\t%1.4f\n', data);

78

fclose(fid);

79 80 81

% Return answers

82

CoV = dimensions(:,2);

83

E = dimensions(:,3);

72

C.6

Global wavelet energy spectrum

The following Matlab code calculates the global wavelet energy spectrum of a slice of each image in a set.

1

% Wavelet analysis script − Richard Riley, 5th Feb 2009

2 3

% GWES = waveline(x real, y real, x pos) calulates the global wavelet

4

% energy spectrum (GWES) of a slice of each selected image (via a dialogue) . The

5

% slice is taken in the y direction (vertical) at position x=x pos where

6

% x pos is in chosen dimensions. x and y real are the real dimensions of

7

% the images in chosen units.

8

function spec = wavline(x real,y real,x pos)

9 10

% Show a dialogue to select the required files.

11

[fname,pname] = uigetfile('*.jpg;*.JPG;*.bmp','Select containing folder','/ home/rriley/Documents/Work/Year 4/Final Project/MatLab/Data/',' multiselect','on')

12

fname = sort(fname);

13 14

% For each image, perform calculations.

15

for k=1:length(fname);

16

f=cell2mat(fname(k))

17

p=([pname f]);

18

rgb image=imread([pname f]);

19 20

% Convert to HSV, to greyscale based on saturation value.

21

conc = im2conc(rgb image,'scale',0);

22

conc show = conc − min(min(conc));

23

conc show = conc show ./ max(max(conc show));

24

figure(1); imshow(conc show); hold on;

25

73

26

% Set the scale range

27

[a,b] = size(conc);

28

a=a−1;

29

scales = 1:a;

30 31

% Define scale factors to give results in real dimensions

32

sfx=double(x real/b);

33

sfy=double(y real/a);

34 35

% Create scale array in real dimensions

36

s real=(1:a).*sfy;

37 38

% Define x pos in terms of pixels

39

x=ceil(x pos/sfx);

40

plot([x,x],[0,a]);

41 42

% Take slice and calculate continuous wavelet transform

43

slice = conc(:,x);

44

figure(2); plot(slice);

45

figure(3);

46

tslice=cwt(slice,scales,'mexh');% ,'plot');

47 48

% Calc GWES

49

for sc=1:a

50

h=tslice(sc,:).ˆ2;

51

spec(sc,2)=(intdump(h,length(h)))*(a*sfy)/(sc*sfy);

52

spec(sc,1)=sc*sfy;

53

end

54 55

% Plot and save to file in selected folder.

56

figure(4); plot(spec(:,1),spec(:,2));

57

fid = fopen([pname f '−wav.dat'], 'w');

58

data=spec';

59

fprintf(fid, '%1.4f\t%1.4f\n', data);

74

60

fclose(fid);

61

end

62

end

75