Modeling of experimental studies of fluid and particle transport in porous media

(1)

Modeling of experimental studies of fluid and

particle transport in porous media

Mikael Härlin Lennermark – Fysikprogrammet Akademin för naturvetenskap och teknik

2009-06-08

Examensarbete, 30 högskolepoäng Handledare: Andreas Oberstedt

(2)

(3)

Abstract

To extract metals and increase the pH value of water around a historical mine waste deposit a series of barrels are used. Polluted water is forced to pass inside these barrels where different filter materials purify the water. This research project is carried out in Sweden by MTM at Örebro University and Bergkraft in Kopparberg, titled “Methods for characterisation and remediation of historical mine waste”. The fluid flow trough the filter materials in the barrels are needed to be understood, in order to improve the extracting process.

In this work a small transparent model filled with sand was made to visualise the fluid flow. In that model coloured water is representing the polluted water. To describe the flow in the transparent model a mathematical model is presented. The theory used in this work is the complex variable method in fluid dynamics together with numerical methods and computer programming. There is a pretty good match between the theoretical and experimental results presented in two dimensions. Continuing work could result in a three dimensional model with different geometries using the same technique.

(4)

(5)

1 Introduction

Historical mining waste sites pollute both surface and ground water with metals and decrease the pH value below acceptable levels. This environmental issue affects the surrounding nature.

At Örebro University the MTM research project “Methods for characterisation and remediation of historical mine waste” in cooperation with the development project Bergskraft Bergslagen is looking for an efficient solution to these problems. The polluted water is forced to pass a series of barrels, which represent filters. Inside the barrels there are several materials, which the metals are extracted onto and the pH value is increased.

The MTM research test site is located at Kopparberg, northern Svealand, Sweden. In Bergslagen alone there are several hundreds of historical mine waste deposits

constantly releasing this acid metal water. This small geographical area in Sweden is not an exception to the rest of the world, which puts this environmental issue in perspective. In Sweden, historical mine waste is one of the most important sources of environmental pollution with trace elements.

In this work a small transparent model was made out of glass to represent the barrels. To simplify the flow pattern the glass model was made rather thin to restrict the problem in nearly two dimensions. Instead of using the inhomogeneous filter material, light coloured sand was used. The experimental method is a picture analysis of

coloured water flowing through the small model. The coloured water is representing the polluted acid metal water.

To better understand the particle transport, which is crucial for the extraction process, a mathematical theory is presented. This theory employs the complex variable method used in fluid dynamics.

After this introduction, theory and experiment are described, followed by the

presentation and discussion of the results obtained in this work. Computer codes and error estimates are given in the Appendix, before I finish with an acknowledgement and a list of references.

(8)

2 Flow Theory

Complex functions

In this chapter a method for solving complicated two dimensional problems in fluid dynamics is presented. The common way of describing a complex variable is,

iy x

z= + , (2.1)

where x and y denote the real and imaginary parts, respectively.

Now w is a complex potential function, with one real and one imaginary function, Ψ

+ Φ

= i

w . (2.2)

If Φ and Ψ are functions of x and y, then w is a function of x and y. Φ is called the velocity potential and Ψ the stream function [1].

The real and imaginary parts of the complex function w, around a line source are given by

and

where r and θ are the polar coordinates, i.e.

(

2 2

)

12

y x r = + and x y arctan = θ .

In two dimensions m is the volume flow rate per unit depth. The depth is referring to the depth of the line source. The fluid is flowing radially from the line source and is therefore seen as a point source in two dimensions. The SI-unit of m is [m2

/s]. For a sink, m is negative [2,3].

The proof of equations (2.3) and (2.4) comes from the Cauchy-Riemann equations for Φ and Ψ, which are necessary for the existence of the function w [4]:

y x ∂ Ψ ∂ = ∂ Φ ∂ (2.5) x y ∂ Ψ ∂ − = ∂ Φ ∂ (2.6) r m ln 2π = Φ _(2.3) θ π 2 m = Ψ , _(2.4)

(9)

Equations (2.3) and (2.4) must be valid for (2.5) and (2.6), which is shown below for the line source mentioned above:

(

)

(

)

(

)

(

)

(

)

(

)

2 2 2 1 2 2 2 1 2 2 2 1 2 2 2 1 2 2 2 1 2 2 2 2 2 1 1 2 1 2 ln 2 ln 2 y x x m x y x y x m y x x y x m y x m x r m x x + ⋅ = ⋅ + ⋅ + ⋅ = =       + ∂ ∂ + ⋅ =       + ∂ ∂ =       ∂ ∂ = ∂ Φ ∂ −

π

(2.7) 2 2 2 2 2 2 2 2 1 2 1 1 1 2 1 1 2 arctan 2 2 y x x m x y x m x x y m x y y x y m x y m y m y y + ⋅ = + ⋅ = ⋅ + ⋅ = =       ∂ ∂ ⋅ + ⋅ =       ∂ ∂ =       ∂ ∂ = ∂ Ψ ∂

π

θ

π

(2.8)

So equations (2.7) and (2.8) meet the first condition (2.5) of the Cauchy-Riemann

equations and the following equations show that the second condition (2.6) is met as

well:

(

)

(

)

(

)

(

)

(

)

(

)

2 2 2 1 2 2 2 1 2 2 2 1 2 2 2 1 2 2 2 1 2 2 2 2 2 1 1 2 1 2 ln 2 ln 2 y x y m y y x y x m y x y y x m y x m y r m y y + ⋅ = ⋅ + ⋅ + ⋅ = =       + ∂ ∂ + ⋅ =       + ∂ ∂ =       ∂ ∂ = ∂ Φ ∂ −

π

(2.9) 2 2 2 2 2 2 2 2 1 1 2 1 1 2 arctan 2 2 y x y m x y x y m x y x x y m x y m x m x x + ⋅ − =       − ⋅ + ⋅ = =       ∂ ∂ ⋅ + ⋅ =       ∂ ∂ =       ∂ ∂ = ∂ Ψ ∂

π

θ

π

(2.10)

Equations (2.3) and (2.4) can be used for an analytical function to describe the flow around a line source or line sink.

Flow around a source

The equation (2.3) describes the potential around a source. Equipotential lines are plotted as contour lines to visualise the potential. The partial derivates of equation (2.3) will describe the velocity around that source, and is recognized as the gradient of Φ [5].

(10)

Using the results of (2.7) and (2.9) the velocity are:

with

and

The equipotential lines around a source (2.3) and the vector plot of the velocity field (2.12) together with (2.13) are plotted below. The vectors are normalized, i.e. transferred into unit vectors. Without the normalization, the difference in magnitude of the vectors is too large to be shown graphically.

Note that the vectors always are perpendicular to the equipotential lines.

Figure 2.1: contour plot of equation (2.3) Figure 2.2: vector plot of equations (2.12) and (2.13)

The contour plot of the stream function (2.4) also intersects at right angles with the contours of equation (2.3) [6].

Figure 2.3: contour plot of equation (2.4) Figure 2.4: surface plot of the potential (2.3)

(

vx vy

)

y x v , _= ,      ∂ Φ ∂ ∂ Φ ∂ = _(2.11) 2 2 2 x y x m x v_x + ⋅ = ∂ Φ ∂ =

π

(2.12) 2 2 2 x y y m y v_y + ⋅ = ∂ Φ ∂ =

π

. (2.13)

(11)

The function (2.3) is proven here to be harmonic, i.e. a solution of the Laplace

equation_{, which is important for the existence of the velocity potential Ф [7]} 0 2 2 2 2 2 = ∂ Φ ∂ + ∂ Φ ∂ ≡ Φ ∇ y x . (2.14)

In the first step, equations (2.7) and (2.9) are used:

(

)

( )

(

)

(

₂ ₂

)

2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 2 2 ` y x x y m y x x x y x m y x x m x x ₊ − = + ⋅ − + ⋅ ⋅ =       + ∂ ∂ = ∂ Φ ∂

π

(2.15)

(

)

( )

(

)

(

₂ ₂

)

2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 2 2 ` y x y x m y x y y y x m y x y m y y ₊ − = + ⋅ − + ⋅ ⋅ =       + ∂ ∂ = ∂ Φ ∂

π

(2.16)

(

) (

)

0 2 ₂ ₂ 2 2 2 2 2 2 2 2 2 2 2 2 2 =         + − + + − = ∂ Φ ∂ + ∂ Φ ∂ = Φ ∇ x y y x x y x y m y x

π

(2.17)

The sum of several harmonic functions is also harmonic, which will be of importance later.

Flow near a wall

The influence of a wall forces the streamlines to bend, so that no particles can cross the wall. The streamlines must then be parallel to the wall. An exact image of the source or sink is used. At half the distance between them a wall is located. This is called the method of images [8].

The two sources need to have equal magnitude of m, here the source is located at (x,y) = (a,b)

(

ln 1 ln 2

)

2 r r m + = Φ

π

, (2.18) where r1 and r2 is

(

) (

)

(

₂ ₂

)

12 1 x a y b r = − + −

(

)

(

)

(

₂ ₂

)

12 2 x a y b r = + + + . (2.19)

Then the velocity components are given by

(

)

(

)

(

) (

)

        ₊ + − = =       ₊ ⋅ ⋅ + − ⋅ ⋅ =       ∂ ∂ ⋅ + ∂ ∂ ⋅ = ∂ Φ ∂ = 2 2 2 1 2 2 1 1 2 2 1 1 2 2 2 1 1 2 2 1 1 2 1 1 2 r a x r a x m r a x r r a x r m x r r x r r m x vx

π

_(2.20)

(12)

and

(

)

(

)

=       ₊ ⋅ ⋅ + − ⋅ ⋅ =       ∂ ∂ ⋅ + ∂ ∂ ⋅ = ∂ Φ ∂ = 2 2 1 1 2 2 1 1 2 2 1 1 2 2 1 1 2 1 1 2 r b y r r b y r m y r r y r r m y v_y π π

(

) (

)

        ₊ + − = ₂ 2 2 1 2 r b y r b y m

π

. (2.21)

The equipotential lines of equation (2.18) with a = 3 and b = 0 and the vector plot of equations (2.20) and (2.21) are shown below. The vectors are normalized. Note that the vectors are parallel to the wall at x = 0.

Figure 2.5: contour plot of (2.18) Figure 2.6: normalized vectors _{(2.20) and (2.21)}

The equation (2.18) is the velocity potential and the stream function is then

(

)

      + + + − − = + = Ψ a x b y a x b y m m arctan arctan 2 2π θ1 θ2 π . (2.22)

To prove that (2.18) can be used, both (2.18) and (2.22) must be valid for the

Cauchy-Riemann equations. The following equations (2.23) to (2.26) are very similar to (2.7)

to (2.10), hence they are presented here in a shorter form:

(

)

_       ₊ + − ⋅ =       ∂ ∂ ⋅ + ∂ ∂ ⋅ ⋅ =       + ∂ ∂ = ∂ Φ ∂ 2 2 2 1 2 2 1 1 2 1 2 1 1 2 ln ln 2 r a x r a x m x r r x r r m r r m x x π π π (2.23)

(13)

(

)

        ₊ + − = =                     + + ∂ ∂ ⋅       + + + +       − − ∂ ∂ ⋅       − − + = =             + + + − − ∂ ∂ =       + ∂ ∂ = ∂ Ψ ∂ 2 2 2 1 2 2 2 1 2 1 1 1 1 2 arctan arctan 2 2 r a x r a x m a x b y y a x b y a x b y y a x b y m a x b y a x b y m y m y y

π

θ

π

(2.24)

So equations (2.23) and (2.24) are identical and then the first condition is met. Next is the proof of the second condition:

(

)

(

) (

)

_       ₊ + − ⋅ =       ∂ ∂ ⋅ + ∂ ∂ ⋅ ⋅ =       + ∂ ∂ = ∂ Φ ∂ 2 2 2 1 2 2 1 1 2 1 2 1 1 2 ln ln 2 r b y r b y m y r r y r r m r r m y y

π

(2.25)

(

)

        ₊ + − − = =                     + + ∂ ∂ ⋅       + + + +       − − ∂ ∂ ⋅       − − + = =             + + + − − ∂ ∂ =       + ∂ ∂ = ∂ Ψ ∂ 2 2 2 1 2 2 2 1 2 1 1 1 1 2 arctan arctan 2 2 r b y r b y m a x b y x a x b y a x b y x a x b y m a x b y a x b y m x m x x

π

θ

π

(2.26)

This implies that there is no problem of adding more points to visualise more walls. Because derivatives of functions with more terms than in equations (2.18) and (2.22) will also meet the conditions of the Cauchy-Riemann equations, this method of adding more functions Φ is sometimes called the superposition principle. The sum will also be solutions of the Laplace equation [9,10].

The following equations imply the statement above, showing that equation (2.18) is harmonic:

(

)( )

(

)

(

)( )

(

)

(

)

(

)

2 2 2 2 1 2 2 3 2 3 1 1 2 2 1 2 2 2 2 4 2 4 1 1 2 2 2 1 2 1 2 2 2 2 1 1 2 m x a x a x x x x r r x a x a m r x a r r x a r r r x a x a m r r r r

π

− − − −    ∂ Φ ∂ ∂Φ ∂ − + =  = _  + _= ∂ ∂  ∂  ∂     ₋   ₊  = _ ⋅ + − −  + ⋅ + + − _=        ₋ ₊  = _ − + − _   (2.27)

(14)

Equation (2.28) is solved similarly as (2.27):

(

)

(

)

2 2 2 2 1 2 2 2 2 4 2 4 1 1 2 2 2 2 2 1 1 2 m y b y b y y y y r r y b y b m r r r r

π

     ∂ Φ ∂ ∂Φ ∂ − + =  = _  + _= ∂ ∂  ∂  ∂     ₋ ₊  = _ − + − _   (2.28)

Adding equations (2.27) and (2.28) leads to

(

)

(

)

(

)

(

)

(

)

(

)

(

)

(

)

(

)

2 2 2 2 2 2 2 2 2 2 4 2 4 2 4 2 4 1 1 2 2 1 1 2 2 2 2 2 2 2 2 4 4 1 2 1 2 2 2 1 2 2 2 4 4 2 2 2 2 1 2 1 2 1 2 1 2 2 2 2 2 1 1 1 1 2 2 2 2 2 2 1 1 1 1 1 1 0 x y x a x a y b y b m r r r r r r r r x a y b x a y b m r r r r r r m m r r r r r r r r

π

∂ Φ ∂ Φ ∇ Φ = + = ∂ ∂  ₋ ₊ ₋ ₊  = _ − + − + − + − _=    ₋ ₊ ₋ ₊ ₊ ₊    = _ + − − _=         =  + − − =  + − − =     ˙ (2.29)

Flow in a 2D model

Similar to equations (2.18) to (2.26) a flow pattern near a right angle corner can be visualised. To simulate a right angle corner, four sources (or sinks) with equal magnitude m have to be used. They are distributed as shown in Fig. 2.7.

Figure 2.7: Source distribution for a right angle corner The velocity potential is then

(

ln 1 ln 2 ln 3 ln 4

)

2 r r r r m + + + = Φ

π

, (2.30) where

(

)

(

)

(

)

(

)

(

)

(

₂ ₂

)

12 2 2 1 2 2 1 b y a x r b y a x r + + + = − + − =

(

) (

)

(

₂ ₂

)

12 3 x a y b r = − + +

(

) (

)

(

₂ ₂

)

12 4 x a y b r = + + − . (2.31)

(15)

The constants a and b are the distances to the walls and so a right angle wall is simulated.

Figure 2.8: contour plot of (2.30) Figure 2.9: The normalized gradient vectors of (2.30)

To simulate a flow inside a two dimensional box a grid of several mirror points are needed. One possible way is to place them out as depicted in Fig. 2.10.

Figure 2.10: The box is shown as a black rectangle.

The corresponding contour plot of the potential Φ with the points distributed as in fig. 2.10 was not symmetric (see figure 2.11).

(16)

Figure 2.11: Asymmetric potential corresponding to Fig 2.10, a = b = 4 cm.

In order to get the potential symmetric around the walls, sources and sinks located at y = -54 and 104 were removed, i.e. the most left and most right points. This resulted in the plot shown in Fig. 2.12.

Figure 2.12: Symmetric potential, a = b = 4 cm. Figure 2.12 shows a contour plot of the potential

        − ⋅ ⋅ = Φ

∑

= = 24 1 24 1 ln ln 1 2 j j i i r r n m

π

, (2.32)

where ri is referring to the sources and rj is referring to the sinks.

Equation (2.32) was solved numerically in Octave and the result is figure 2.12. The 48 functions of r(x,y) are showed in Appendix A.1, the source code for the Octave [11] programming is given in Appendix A.2.

The fluid flows in a porous medium and its porosity n has to be taken into account. It is defined as the total volume of the pores divided by the total (bulk) volume [12]. The strength of the source m is related to the real three-dimensional volume flow q as

(17)

d t V d q m= = ⋅1 , _(2.33)

where d corresponds to the thickness of the volume that is considered. The volume available for the flow is restricted to the volume of the pores, which is why the porosity n is an important factor. By definition n is a number in the interval 0 < n < 1. The velocity of the fluid is the partial derivative of the potential described by equation (2.32),

and

In the above equations, x – CX and y – CY are results of the derivatives of r. The constants CX and CY can be seen in Appendix A.1 in the functions r, expressed in terms of H, L, a and b, where H is the height and L is the length of the model volume.

Figure 2.13: Equipotential lines and the streamlines are showed in the area of the model, a = b = 0.1 cm

The streamlines in figure 2.13 are by definition, tangential to the velocity of the fluid [13], (compare figures 2.13 and 2.14).

        ₋ − − ⋅ ⋅ = ∂ Φ ∂ =

∑

= = 24 1 2 24 1 2 1 2 j _j j j i _i i i x r CX x r CX x n m x v

π

(2.34)         ₋ − − ⋅ ⋅ = ∂ Φ ∂ =

∑

= = 24 1 2 24 1 2 1 2 j _j j j i _i i i y r CY y r CY y n m y v

π

. (2.35)

(18)

Figure 2.14: Normalized velocity of the fluid

Figure 2.14 shows the vector plot using equations (2.34) and (2.35). Appendix A.3 shows the Octave source code, where equations (3.34) and (2.35) are expressed as an Octave function.

Path lines

By using the simple relationds= v⋅dt, or expressed for the numerical approach as

model

s=

∑

v⋅ ∆t

,

(2.36) the path lines may be visualised. Here s is the distance of the fluid particle, v is its velocity and t is the time. In this stationary case path lines and streamlines coincide, compare Figures 2.15 and 2.13. The source code is seen in part 4 of Appendix A.3.

(19)

Figure 2.16 shows locations of the fluid particles after different times t. The strength of source m has for the experimental study typical value of m = 8.8 · 10-3 cm2s-1 and the distance of the sources/sinks to the walls is a=b=0.5cm. The Octave source code used to plot these graphs is seen in Appendix A.4.

t = 0,5 h

t = 1 h

t = 8 h

t = 9 h

(20)

3 Experiments

The objectives of the MTM research project are to extract metals from the polluted water and increase its pH value. These impurities along with acidity are an

environmental issue around historical mining waste sites. A series of filters as seen in figure 3.1 is one probable solution to this problem.

In this work a method of visualising the fluid flow inside the filter barrels of figure 3.1 are presented.

Figure 3.1: Photograph of the site near Kopparberg.

Transparent model

To visualise a flow through wet sand a transparent model was made out of glass. The length was 50 cm and the height 20 cm. To force the fluid in mainly two dimensions the box was only 5 cm wide.

Sand in the grain size interval of 0.1 mm to 0.3 mm and another type with grains in the interval of 0.4 mm to 0.8 mm was used in different experiments.

Figure 3.2: Glass box for the experimental simulation of the particle flow in sand

(21)

Upper reservoir

Lower reservoir

Constant water level Tap to regulate the flow

Flexible tube

Setup

A constant flow throughout the experiment was needed, hence two reservoirs were used. The upper reservoir was only opened at the bottom, as seen in figure 3.4, whereas the water level of the bottom reservoir was stabilised. The lower reservoir was open at the top and at the bottom a flexible plastic tube was attached. It supplies the box with water, whose flow could be regulated by a tap. Since the pressure due to the water volume in the upper vessel is compensated by the air pressure inside, the water level in the lower reservoir is only exposed to atmospheric pressure. When the water level decreased below the opening of the upper reservoir, water was refilled from there until the opening is filled again with water. This way the water level varies only within a few millimetres, resulting in a nearly constant flow. Hence, the choice of using this setup with two coupled reservoirs together with a tap provides an excellent stabilisation of the flow rate.

The total setup is shown in figure 3.3, while details are depicted in Figure 3.4.

Throughout the experiments the inclination of the box was set at 10° to represent the inclination at the MTM test site. In order to detect the water flow a digital camera on a stand and a clock were used. To visualise the movement of the fluid through the wet sand a colour was inserted in the water. The colour used in the experiments is

completely soluble in water. Figure 3.5 shows how easily the fluid could be observed.

Figure 3.5: Coloured water flow through wet sand

Figure 3.3: Setup overview

(22)

The pictures were analysed by means of the computer programme Microsoft Paint™ [14] which has a proper coordinate system. The path of the colour was then able to be compared with the numerical results. To calculate the volume flow rate the coordinate system was used to measure the increasing volume in the water cup seen at the left end in figure 3.5. The experimental values are summarised in Appendix A.5.

The porosity was calculated by measuring the volume of the sand and the mass of the water that completely filled the pores, i.e. making the sand completely saturated. In Appendix A.6 the data and the calculation of the porosity are shown.

4 Results and discussion

On the following pages (Figure 4.1 to 4.5) pictures taken at different times are compared with the numerical results. The fluid path s is plotted together with its uncertainty of 17%. The uncertainty estimate is discussed later.

The following pictures exhibit some difference between experimental and numerical results that may be explained with the not entirely constant flow rate. The method of keeping the flow rate constant could be improved. The flow rate per unit depth m in the numerical results is based upon the average of the experimental flow rate q and the average is calculated over some hour or hours. In this long time span some fluctuations in the flow rate q are possible.

When the graph is plotted on top of the picture the edge of the coloured and non-coloured water was sometimes difficult to see. Hence the changing in contrast and colour intensity was necessary, which was performed with Adobe Photoshop™ [15].

(23)

Figure 4.1 t = 1h n = 0,54 _{± 5%} q = 0,089cm 3_/s ± 12% m = 0,018 cm2/s ± 12%

(24)

Figure 4.2 t = 2,5 h n = 0,54 ± 5% q = 0,089cm3/s ± 12% m = 0,018 cm2/s ± 12%

(25)

Figure 4.3 t = 4 h n = 0,54 _{± 5%} q = 0,089cm 3_/s ± 12% m = 0,018 cm2/s ± 12%

(26)

Figure 4.4 t = 20 minutes n = 0,50 ± 4% q = 0,28cm3/s ± 13% m = 0,056 cm2/s ± 13%

(27)

Figure 4.5 t = 45 minutes n = 0,50 ± 4% q = 0,28cm3/s ± 13% m = 0,056 cm2/s ± 13%

(28)

Error estimate

The flow was calculated from analysing the increasing water level in the water cup seen at the bottom left in figures 4.1 to 4.5. To estimate the error the standard deviation was used. The standard deviation was 12% from the average in the experiment concerning figures 4.1 to 4.3 and 13% in the other one, figures 4.4 and 4.5, see Appendix A.5.

The porosity was found by determine the bulk density of dry sand and completely saturated sand. To measure the porosity the mass of the water in the wet sand was also known (see Appendix A.6). The standard deviation of the porosity was 5% from the average in the sand with grain size 0.1 mm to 0.3 mm. In the sand with grain size 0.4 mm to 0.8 mm the standard deviation was 4%.

Both the error concerning the flow and the porosity could of course decrease with more measurements and more accurate equipment, which however was not possible within the time frame of this work.

So the uncertainties of the flow and the porosity were easily found, but there are other errors involved. To get the numerical results the distribution of the sand grains is not considered, only the porosity. The sand tended to distribute in layers. Therefore it is important to mix the sand properly. If the sand grains were perfectly spherical, mixing would not be necessary. Experiments showed that, if the box was filled without mixing, the water followed the sand layers. And the experiment would give

completely other results than the numerical ones. These layers are seen in figure 4.6, where the water is flowing through dry sand. Here the layers are nearly horizontal.

Figure 4.6: water flowing in sand layers

Another error involved is the fact that the flow through the sand is only visible from one side. What the fluid path looks like further inside of the model is not seen. This problem is minimized by using a rather thin box.

Using the fractional uncertainties of the flow and the porosity the fractional uncertainty of the fluid path is

s m n q n

δ

≤ + = + (4.1)

The volume flow rate q and the flow rate per unit depth m is proportional as m=q/d, where d is the width of the model and so the relative errors of m and q are equal.

(29)

The fractional uncertainty of the path s in the experiment concerning figures 4.1 to 4.3 and 4.4 to 4.5 are both 17%. This is shown by error bars in these figures.

About this method

With this method, this work has proven to be powerful when solving different problem regarding flows. In order to study another fluid path than in the transparent model one could change the location of the sources, sinks and mirror points. The fluid path is even more accurate with more mirror points. And without computer

programming none of these results are possible to obtain.

Ideas of more work

There is a missing part of the flow in the numerical results: the excess flow at the surface of the sand. In the pictures, that part is visible as a tail near the sand surface. To correct this, which would be possible, one way is to decrease m somewhat and let that excess part of m flow at the surface. This is resulting in a new source that would become a line source in two dimensions, starting from the real point source and increase its length as time goes on. This line source would add extra fluid from the surface, creating a tail. Moreover, an adjustment of the porosity of the sand to the experimental observations might lead to an improved agreement between theoretical description and experiment.

The aim of this work was to describe a two dimensional flow. However, it would also be possible in three dimensions. The velocity potential and the stream function follow contour lines in two dimensions. In three dimensions the velocity potential is seen as equipotential surfaces perpendicular to the stream functions surfaces [16].

(30)

(31)

Appendix

A.1 The functions r from equation (2.32)

Sources (index i)

(

)

(

)

(

₂ ₂

)

12 1 x a y 2H b r = − + − −

(

) (

)

(

₂ ₂

)

12 2 x a y 2H b r = − + − +

(

) (

)

(

₂ ₂

)

12 3 x a y 2H b r = + + − −

(

) (

)

(

₂ ₂

)

12 4 x a y 2H b r = + + − +

(

) (

)

(

₂ ₂

)

12 5 x a y H b r = − + − −

(

) (

)

(

₂ ₂

)

12 6 x a y H b r = − + − +

(

) (

)

(

₂ ₂

)

12 7 x a y H b r = + + − −

(

) (

)

(

₂ ₂

)

12 8 x a y H b r = + + − +

(

) (

)

(

₂ ₂

)

12 9 x a y H b r = − + + −

(

) (

)

(

₂ ₂

)

12 10 x a y H b r = − + + +

(

) (

)

(

₂ ₂

)

12 11 x a y H b r = + + + −

(

) (

)

(

₂ ₂

)

12 12 x a y H b r = + + + +

(

) (

)

(

₂ ₂

)

12 13 x a y 2H b r = − + + −

(

)

(

)

(

₂ ₂

)

12 14 x a y 2H b r = − + + +

(

) (

)

(

₂ ₂

)

12 15 x a y 2H b r = + + + −

(

) (

)

(

₂ ₂

)

12 16 x a y 2H b r = + + + +

(

)

(

)

(

₂ ₂

)

12 17 x 2L a y 2H b r = − + + − −

(

)

(

)

(

₂ ₂

)

12 18 x 2L a y 2H b r = − + + − +

(

) (

)

(

₂ ₂

)

12 19 x 2L a y H b r = − + + − −

(

) (

)

(

₂ ₂

)

12 20 x 2L a y H b r = − + + − +

(

) (

)

(

₂ ₂

)

12 21 x 2L a y H b r = − + + + −

(

) (

)

(

₂ ₂

)

12 22 x 2L a y H b r = − + + + +

(

)

(

)

(

₂ ₂

)

12 23 x 2L a y 2H b r = − + + + −

(

) (

)

(

₂ ₂

)

12 24 x 2L a y 2H b r = − + + + + Sinks (index j)

(

) (

)

(

₂ ₂

)

12 1 x L a y 2H b r = − − + − −

(

) (

)

(

₂ ₂

)

12 2 x L a y 2H b r = − − + − +

(

) (

)

(

)

12 3 x L a y 2H b r = − + + − −

(

) (

)

(

)

12 4 x L a y 2H b r = − + + − +

(

)

(

)

(

₂ ₂

)

12 5 x L a y H b r = − − + − −

(

)

(

)

(

₂ ₂

)

12 6 x L a y H b r = − − + − +

(

)

(

)

(

₂ ₂

)

12 7 x L a y H b r = − + + − −

(

)

(

)

(

₂ ₂

)

12 8 x L a y H b r = − + + − +

(

) (

)

(

₂ ₂

)

12 9 x L a y H b r = − − + + −

(

)

(

)

(

₂ ₂

)

12 10 x L a y H b r = − − + + +

(

) (

)

(

₂ ₂

)

12 11 x L a y H b r = − + + + −

(

) (

)

(

₂ ₂

)

12 12 x L a y H b r = − + + + +

(

)

(

)

(

₂ ₂

)

12 13 x L a y 2H b r = − − + + −

(

) (

)

(

₂ ₂

)

12 14 x L a y 2H b r = − − + + +

(

)

(

)

(

₂ ₂

)

12 15 x L a y 2H b r = − + + + −

(

)

(

)

(

₂ ₂

)

12 16 x L a y 2H b r = − + + + +

(

)

(

)

(

₂ ₂

)

12 17 x L a y 2H b r = + + + − −

(

)

(

)

(

₂ ₂

)

12 18 x L a y 2H b r = + + + − +

(

)

(

)

(

₂ ₂

)

12 19 x L a y H b r = + + + − −

(

)

(

)

(

₂ ₂

)

12 20 x L a y H b r = + + + − +

(

)

(

)

(

₂ ₂

)

12 21 x L a y H b r = + + + + −

(

)

(

)

(

₂ ₂

)

12 22 x L a y H b r = + + + + +

(

) (

)

(

₂ ₂

)

12 23 x L a y 2H b r = + + + + −

(

) (

)

(

₂ ₂

)

12 24 x L a y 2H b r = + + + + +

(32)

A.2 potential.m

Part 1 Part 2

% octave script to calculate and plot equipotential lines % one source and one sink in a 2D-model

% april 2009

% Mikael Härlin Lennermark

nr = 200; % creating vectors and a grid x = linspace(-60,110,nr);

y = linspace(-60,60,nr); [X,Y] = meshgrid(x,y);

L = 50; % lenght and height of the box [cm] H = 17; Nj = 4; Ni = 4; % assigning 16 points for j = 1:Nj for i = 1:Ni xpoint(i, j) = i * L; ypoint(i, j) = j * 2 * H; endfor endfor

% to locate the real source and sink at the desired place xpoint = xpoint .- 2 * L;

ypoint = ypoint .- 5 * H;

m = 1; %volume flow rate per unit depth [cm^2 s^-1] a = 4; % a and b is the distance source to walls[cm] b = 4;

Nk = (4 * Ni * Nj); xterm = zeros(nr,nr,Nk); yterm = zeros(nr,nr,Nk);

% giving coordinates to the source, sink and all mirror % points k = 1; while (k < Nk for j = 1:Nj for i = 1:Ni xterm(:, :, k) = X .- xpoint(i, j) .- a; xterm(:, :, k+1) = X .- xpoint(i, j) .+ a; xterm(:, :, k+2) = X .- xpoint(i, j) .- a; xterm(:, :, k+3) = X .- xpoint(i, j) .+ a; yterm(:, :, k) = Y .- ypoint(i, j) .+ b; yterm(:, :, k+1) = Y .- ypoint(i, j) .+ b; yterm(:, :, k+2) = Y .- ypoint(i, j) .- b; yterm(:, :, k+3) = Y .- ypoint(i, j) .- b; k = k + 4; endfor endfor endwhile

% r() is the radial distance to the 64 (=Nk) points r = zeros(nr,nr,Nk); for k = 1 : Nk r(:, :, k) = ((xterm(:, :, k)).^2 .+ (yterm(:, :, k)).^2).^(1/2); endfor % the potential potential = zeros(nr); k = 1; while (k < Nk)

for k1 = k : (k+3) % sinks in this interval

if (k1 != 2 & k1 != 4 & k1 != 18 & k1 != 20 & k1 != 34 & k1 != 36 & k1 != 50 & k1 != 52)

potential = potential .- log(r(:, :, k1)); endif

endfor

for k2 = (k+4) : (k+7) % sources in this interval

if (k2 != 13 & k2 != 15 & k2 != 29 & k2 != 31 & k2 != 45 & k2 != 47 & k2 != 61 & k2 != 63)

potential = potential .+ log(r(:, :, k2)); endif endfor k = k + 8; endwhile

potential = (m/(2*pi)) * potential;

%plotting the equipotential lines in black ('k') contour(X, Y, potential, 50,'k')

xlabel('x [cm]') ylabel('y [cm]')

(33)

A.3 particles_path.m

Part 1 Part 2

% octave script to calculate coordinates for particle % path in a flow and write it in a text file

% one source and one sink in a 2D-model % april 2009

% Mikael Härlin Lennermark

L = 50; % lenght and height of the box [cm] H = 17;

a = 0.5; b = 0.5;

% creating an array of angels for the particles start % points

particle_nr = 40; deg = 0;

for i = 1:particle_nr degrees(i) = deg; deg = deg + 360/particle_nr;

endfor

teta = degrees * (2*pi/360); % in radians % creating an array of cartesian coordinates for the %particles start points

distance_from_source = 0.1; for i = 1:particle_nr

x_start(i) = a + distance_from_source * cos(teta(i)); y_start(i) = H - b - distance_from_source * sin(teta(i)); endfor

% a function to find the velocity, the partial derivative % of the potential

function [V_X V_Y] = velocity (X, Y)

L = 50; % lenght and height of the box [cm] H = 17;

% the distance source to walls or sink to walls [cm] a = 0.5;

b = 0.5; Nj = 4; Ni = 4;

% assigning 16 points (source and sink + 14 mirror %points) for j = 1:Nj for i = 1:Ni xpoint(i, j) = i * L; ypoint(i, j) = j * 2 * H; endfor endfor

% to locate the real source and sink at the desired place xpoint = xpoint .- 2 * L;

ypoint = ypoint .- 5 * H;

% volume flow rate per unit depth [cm^2 s^-1] m = 0.055845;

% porosity of the sand n = 0.49871; Nk = (4 * Ni * Nj);

xterm = linspace(0,0,Nk); yterm = linspace(0,0,Nk);

%coordinates to the source, sink and all mirror points k = 1; while (k < Nk) for j = 1:Nj for i = 1:Ni xterm(k) = X - xpoint(i, j) - a; xterm(k+1) = X - xpoint(i, j) + a; xterm(k+2) = X - xpoint(i, j) - a; xterm(k+3) = X - xpoint(i, j) + a; yterm(k) = Y - ypoint(i, j) + b; yterm(k+1) = Y - ypoint(i, j) + b; yterm(k+2) = Y - ypoint(i, j) - b; yterm(k+3) = Y - ypoint(i, j) - b; k = k + 4; endfor endfor endwhile

% r() is the radial distance to the 64 (=Nk) points r = linspace(0,0,Nk); for k = 1 : Nk r(k) = ((xterm(k))^2 + (yterm(k))^2)^(1/2); endfor

V_X = 0; % The velocity components V_Y = 0;

(34)

Part 3 Part 4 k = 1;

while (k < Nk)

for k1 = k : (k+3) % sinks in this interval if (k1 != 2 & k1 != 4 & k1 != 18 & k1 != 20 & k1 != 34 & k1 != 36 & k1 != 50 & k1 != 52)

V_X = V_X - (xterm(k1) / (r(k1)^2)); V_Y = V_Y - (yterm(k1) / (r(k1)^2)); endif

endfor

for k2 = (k+4) : (k+7) % sources in this interval

if (k2 != 13 & k2 != 15 & k2 != 29 & k2 != 31 & k2 != 45 & k2 != 47 & k2 != 61 & k2 != 63)

V_X = V_X + (xterm(k2) / (r(k2)^2)); V_Y = V_Y + (yterm(k2) / (r(k2)^2)); endif endfor k = k + 8; endwhile V_X = m/(2*pi*n) * V_X; V_Y = m/(2*pi*n) * V_Y; endfunction

% open a text file, FILE is a pointer FILE = fopen('coordinates2.dat', 'w');

% a loop to print the path of the particles in the text file

for i = 1:particle_nr

t = 0.0; % t is time for each particle x(i) = x_start(i);

y(i) = y_start(i); while (x(i) < L-0.6)

dt = 50 * sin((pi*x(i))/50); % sin() is used to % get short dt at the beginning and at the

% end of the interval

t = t + dt;

fprintf(FILE, '%f\t%f\t%f\n', t, x(i), y(i));

[v_x v_y] = velocity(x(i), y(i)); % calling function

x(i) = x(i) + dt * v_x; % assigning next coordinate y(i) = y(i) + dt * v_y;

endwhile fprintf(FILE, '%s\n', ' '); endfor fclose(FILE);

(35)

A.4 plot_particles.m

% octave script to plot fluid particles after time t % april 2009

% Mikael Härlin Lennermark % set a time [s]

time = 1200;

% tolerance is the interval (time - tolerance) < time < (time + tolerance) tolerance = 15; % in seconds

% error estimation

fractional_uncertainty = 0.17; % open the file with the coordinates FILE = fopen('coordinates2.dat' , 'r'); t = fscanf(FILE, '%f', 1); x = fscanf(FILE, '%f', 1); y = fscanf(FILE, '%f', 1); x_error = x * fractional_uncertainty; y_error = y * fractional_uncertainty;

x = 0 - x; % to get the particles coming from right to left similar with the pictures number_of_plots = 0;

% loop until end of file is found and plot the positions while (! feof(FILE) )

if (t > (time - tolerance) & t < (time + tolerance) ) plot(x, y)

errorbar( x, y, x_error, x_error, y_error, y_error); hold on number_of_plots = number_of_plots + 1; endif t = fscanf(FILE, '%f', 1); x = fscanf(FILE, '%f', 1); y = fscanf(FILE, '%f', 1); x_error = x * fractional_uncertainty; y_error = y * fractional_uncertainty; x = 0 - x; endwhile number_of_plots fclose(FILE);

xlabel('x [cm]') % Labels to the axis

ylabel('y [cm]')

axis ([ -50, 0, 0, 17.5]) % intervals of the axis

axis('equal'); % forcing the relative unit lenght of the axis to be equal hold off

(36)

A.5 Flow Calculation

To measure the flow rate q the increasing water height in the water cup was measured. This was done with several pictures in Microsoft Paint™ with has a proper coordinate system in the units of pixels. This is converted to centimetres. The length of the model is 50.8 cm (*). The water cup is cylindrical with a radius of 5.55 cm (**).

Experiment with flow in sand grain size 0.1 – 0.3 mm (Figure 4.1 to 4.3)

length of model conversion factor * water height water height volume ** time dV Dt flow rate q

[pixels] [pixels/cm] [pixels] [cm] [cm3] [s] [cm3] [s] [cm3s-1]

1470,9 28,96 25 0,86 83,5 3601 1462,7 28,79 40 1,39 134,4 4200 50,9 599,0 0,085 1477,4 29,08 54 1,86 179,7 4808 45,2 608,0 0,074 1470,8 28,95 70 2,42 234,0 5401 54,3 593,0 0,092 1471,5 28,97 100 3,45 334,1 6604 100,1 1203,0 0,083 1471,8 28,97 171 5,90 571,2 3000 1483,2 29,20 119 4,08 394,4 1200 176,8 1800 0,098 1481,4 29,16 100 3,43 331,8 600 62,6 600 0,104 Average [cm3s-1] std. dev [cm3s-1] 0,089 0,011 fractional uncertainty 12% Experiment with flow in sand grain size 0.4 – 0.8 mm (Figure 4.4 to 4.5)

length of model conversion factor * water height water height volume ** time dV Dt flow rate q

[pixels] [pixels/cm] [pixels] [cm] [cm3] [s] [cm3] [s] [cm3s-1]

1502,4 29,57 96 3,25 314,1 0 1505,6 29,64 114 3,85 372,2 196 58,1 196 0,296 1504,5 29,62 124 4,19 405,2 301 32,9 105 0,314 1504,0 29,61 148 5,00 483,7 601 78,6 300 0,262 1504,4 29,61 172 5,81 562,0 903 78,3 302 0,259 1502,4 29,57 34 1,15 111,2 1200 297 1502,4 29,57 57 1,93 186,5 1502 75,3 302 0,249 1504,7 29,62 89 3,00 290,8 1804 104,3 302 0,345 1503,7 29,60 158 5,34 516,5 2698 225,8 894 0,253 1499,4 29,52 62 2,10 203,3 3602 904 1499,6 29,52 132 4,47 432,7 4500 229,4 898 0,256 Average [cm3s-1] std. dev [cm3s-1] 0,279 0,035 fractional uncertainty 13%

(37)

A.6 Porosity calculation

The porosity n is the volume of the pores divided by the total volume. Here the volume of the pores is measured by measuring the maximum mass of water contained in the pores. The mass was found using a ±1g weighing scale and the volume using a ±1ml measuring cup. First the average density of the dry sand is needed for the porosity calculation.

Grain size 0,1-0,3 mm Grain size 0,4-0,8 mm

mdry [g] Vdry [ml] ρdry [kg/m3] mdry [g] Vdry [ml] ρdry [kg/m3] 16 12 1,33E+03 28 19 1,47E+03 27 20 1,35E+03 17 11 1,55E+03 40 30 1,33E+03 24 17 1,41E+03 49 37 1,32E+03 57 38 1,50E+03 63 47 1,34E+03 100 65 1,54E+03 90 68 1,32E+03 136 89 1,53E+03 117 87 1,34E+03 42 29 1,45E+03 130 95 1,37E+03 69 47 1,47E+03 135 100 1,35E+03 89 60 1,48E+03 81 62 1,31E+03 19 12 1,58E+03 41 36 1,14E+03 34 21 1,62E+03 54 35 1,54E+03 62 41 1,51E+03 25 17 1,47E+03 101 68 1,49E+03 21 15 1,40E+03 126 81 1,56E+03 24 19 1,26E+03 152 98 1,55E+03 20 14 1,43E+03 23 16 1,44E+03 Average Std 18 12 1,50E+03 * dev 26 18 1,44E+03 [kg/m3] [kg/m3] 1,51E+03 54 Average * [kg/m3 ] Std dev [kg/m3 ] 1,37E+03 91 Porosity calculation

(

)

₍

₎

tot tot w dry tot tot tot tot w dry tot tot tot w w tot w m V m m m m m m V V n

ρ

1 ⋅ ⋅ − = − = = = where n is the porosity,

Vw The volume of the water inside the pores,

Vtot total volume (dry sand + water),

mw mass of the water, ρw density of the water,

mtot total mass, ρtot total density,

mdry mass of the dry sand and

ρdry * the density of the dry sand e.g. the average from the table above.

Vtot=Vdry This assumption is made. Vdry is the volume of the dry sand. Next is the calculation of the porosity.

(38)

Grain size 0,1-0,3 mm Grain size 0,4-0,8 mm mtot [g] Vtot [ml] ρtot [kg/m3] n [-] mtot [g] Vtot [ml] ρtot [kg/m3] n [-] 45 24 1,88E+03 0,51 30 15 2,00E+03 0,49 94 51 1,84E+03 0,47 50 25 2,00E+03 0,49 136 71 1,92E+03 0,55 78 39 2,00E+03 0,49 149 78 1,91E+03 0,54 97 49 1,98E+03 0,47 167 87 1,92E+03 0,55 140 69 2,03E+03 0,52 175 91 1,92E+03 0,55 175 86 2,04E+03 0,52 184 96 1,92E+03 0,55 187 93 2,01E+03 0,50 38 20 1,90E+03 0,53 192 95 2,02E+03 0,51 63 33 1,91E+03 0,54 55 28 1,96E+03 0,45 68 35 1,94E+03 0,57 86 42 2,05E+03 0,53 89 47 1,89E+03 0,53 105 52 2,02E+03 0,51 106 55 1,93E+03 0,56 115 57 2,02E+03 0,50 121 63 1,92E+03 0,55 128 64 2,00E+03 0,49 150 74 2,03E+03 0,51

Average Std dev 164 81 2,03E+03 0,51

[-] [-] 180 89 2,02E+03 0,51 0,54 0,03 187 93 2,01E+03 0,50

_uncertaintyfractional Average

[-] Std dev [-] 5% 0,50 0,02 fractional uncertainty 4%

(39)

Acknowledgments

First I would like to thank my excellent supervisor Professor Andreas Oberstedt at Örebro University, who has supported me and all the time believed in me. His enthusiasm and devotion to physics has been very important to me.

I am very much grateful to MTM and Associate Professor Stefan Karlsson at Örebro University who has given me the opportunity to work with this project. And thanks to the rest of the team involved in the MTM project: Mattias Bäckström and Lotta Sartz. The powerful tool of numerical methods was brought to my attention by Professor Peter Johansson at Örebro University. His useful methods in Octave have been of great importance to complete this work.

I also want to thank another teacher at Örebro University, lecturer Fredrik Wallinder, who together with Andreas Oberstedt and Peter Johansson has been an inspiration during my physics studies.

My excellent colleagues throughout the years at Örebro University, Robert Billnert, Johan Karlsson and Johannes Morfeldt have been very fun studying with. Our support to each other has been a key to succeed our Physics studies.

Thanks to all other MTM workers and scientific staff at Bilbergska huset, Örebro University, who I spent some time with during coffee breaks and lunchtime. A special thanks to my parents who always have encouraged me.

Last but not least, my gratitude goes to my wonderful family, my wife Lisa and our son Algot, who has been a very important support to me.

(40)

(41)

References

[1] P. K. Kundu and I. M. Cohen, Fluid Mechanics, Elsevier Academic Press, (2004), p. 157-158

[2] E. Boeker and R. Van Grondelle, Environmental Physics, John Wiley & Sons Ltd, (1999), p. 229

[3] P. K. Kundu and I. M. Cohen, Fluid Mechanics, Elsevier Academic Press, (2004), p. 161

[4] E. Boeker and R. Van Grondelle, Environmental Physics, John Wiley & Sons Ltd, England, (1999), p. 224

[6] G. Backstrom, Fluid Dynamics by Finite Element Analysis, Studentlitteratur, (1999), p. 25

[9] P. K. Kundu and I. M. Cohen, Fluid Mechanics, Elsevier Academic Press, (2004), p. 159

[11] J.W. Eaton and others, GNU Octave, version 3.0.3 (2008)

[13] J. F. Douglas et al., Fluid Mechanics, Longman Scientific & Technical, (1995), p. 88

[14] Microsoft Corporation, Microsoft Paint, version 5.1 (2007) [15] Thomas Knoll et al., Adobe Photoshop, version 8.0.1 (2003)

Modeling of experimental studies of fluid and particle transport in porous media