Numerical and experimental studies of annular flume flow

ELSEVIER

Numerical and experimental flume flow

studies of annular

P. W. James, T. E. R. Jones and D. M. Stewart School

of Mathematics

and Statistics,

Unioersity

of Plymouth,

Plymouth,

UK

Annular flume studies of cohesive sediment beds are routinely used to investigate the erosion characteristics of the bed and the transport of suspended sediment. In order to interpret the results of such studies, an accurate description of the f7ow field and bed stress distribution is required. In this paper the fzow field and bed stress distribution are calculated for a miniature annular flume using the computational fluid dynamics (CFD) package Harwell-FLOW3D. Particular attention is paid to the effect of singularities at the edges of the annular lid. The predictions are compared with data obtained from a miniature annular flume mounted on a Carri-med rheometer. By modifiiing the package, calculations are also made for inelastic non-Newtonian fluids, and the results are compared with data obtained with a Weissenberg rheogoniometer and flush mounted hot-film shear stress probes.

Keywords: annular flume, cohesive sediment, flow calculation, non-Newtonian fluids, shear stress measurement

1. Introduction For some time it has been recognized that annular flumes provide a useful means of measuring certain properties of cohesive sediment beds, notably their erosion characteristics.’ The precise mechanisms that govern the erosion and transport of a cohesive sediment bed are still not completely understood because of the complex chemical and physical processes that are involved. However flume studies have been used to clarify some of the dominant mechanisms, and it is accepted that the shear stress at the interface of the cohesive sediment bed and the overlying fluid is an important controlling parameter. There appears to be a critical shear stress above which significant erosion of the sediment bed takes place. Straight flumes have the disadvantage that the suspended sediment has to be recirculated, and the pumps and recirculation loops required can significantly alter its structure. Annular flumes, in which the flow is driven by a rotating lid, overcome this problem and so have become popular in recent years both in the laboratory and in the field.2 However the flow in an annular flume is not unidirectional and so an accurate description of the flow field is required before detailed interpretation of experimental results can take place. In

Address reprint requests to Dr. James at the School of Mathematics and Statistics, University of Plymouth, Drake Circus, Plymouth PL4 EAA, UK. Received

12 August

1993; accepted

2 February

1994

Appl. Math. Modelling 1996, Vol. 20, March 0 1996 by Elsevier Science Inc. 655 Avenue of the Americas, New York, NY 10010

particular the bed stress distribution must be calculated accurately to determine the erosion properties of the sediment bed. In previous work,3 numerical predictions of the flow field and surface stress distribution in the Hydraulics Research (Wallingford) annular flume4 have been compared with experimental results. The numerical predictions were made with Harwell-FLOW3D, using both a k-c and a mixing length turbulence model, the velocity field was measured using laser Doppler anemometry, and surface stresses were measured with hot-film flush-mounted sensors. Very good agreement was obtained between predicted and measured velocity profiles, and good agreement between measured and predicted surface stress profiles was also found (the stresses on the curved annular walls being in closer agreement than the bed stresses). In this paper the bed stress and lid torque are measured and calculated for a miniature annular flume, the miniflume (Figure I), in which it is planned to carry out studies of the erosion and transport of cohesive sediments. James and Jones’ have shown that it is possible to obtain good agreement between theory and experiment for the flow of an almost constant viscosity (0.94 Pa . s) silicone oil and a shear-thinning elastic liquid. Here we present a more detailed numerical investigation of the effect of the singularities at the edges of the lid, and it is shown that the results obtained previously are sufficiently accurate for comparison with experimental data. New experimental results are also presented for a miniature annular flume in which the torque on the fixed flume bed can be measured. Comparison of these results with predictions from Harwell-FLOW3D confirms that the code, modified for use

0307-904X/96/$15.00 SSDI 0307-904X(95)00149-2

Annular

flume

flow:

P. W. James

et al. Here 77= v(q) is the shear-dependent where y is the shear rate, defined by

ROTATING

_-_

fluid

viscosity,

au 2

LID

+z (1 +(;+;)2+(::f)2}*‘2 (6)

i={2(32+2(32+2(32

,e

r-a

I

The boundary

Figure 1. A sketch of the miniflume tem (a=6.4cm,

b=2.95

and the coordinate cm, a=O.l48cm, p=O.l8cm).

conditions

for the above equations

n=~,=w=O at r=a, r= b, Olz
with shear-dependent viscosity, provides description of the bed stress distribution.

2. Mathematical

on z=h,

a very accurate

The problem considered is that of laminar flow in an annulus of depth h driven by a moving lid. With respect to a cylindrical polar coordinate system (I, 0, z), the annulus is defined by a 5 r IS b, 0 I z 5 h (Figure I). The rotating annulus lid occupies the region a + cx I r I b - /3, z = h, where (Y and /3 are the gaps between the edges of the rotating lid and the fixed side walls at r = a and r = b, respectively. The fluid is assumed to be incompressible and the flow to be steady and axisymmetric. The governing equations are then

3. Zero Reynolds

I

number,

Newtonian

flow

In order to assess the influence of different conditions in the gaps between the rotating lid side walls, the problem formulated in Section 2 fied by restricting it to constant viscosity and low number. In this type of problem it is appropriate “rotational” Reynolds number, Re, defined by

dw (1)

r --+az=O dr

(7)

w=rLl a+n
where R is the angular speed of the rotating lid. In the gaps, a < r < a + (Y, b - /3 < r < b, there are free surfaces but, for simplicity, we assume that these surfaces are horizontal (i.e., along the line z = h). The precise boundary conditions to be applied in the gaps are discussed in the next section.

formulation

1 a( ru)

.nb2p Re = P

au

uv

p uar+-+wz

r

i

au

780

JP

r

dr

G(rzTo) +i(TzOl ) :12 (3) (4)

In these equations, (u, v, w) is the velocity vector, p the pressure, p the density, and rrZ etc. are the components of the shear stress tensor and are defined by

U

aw\

au 7 20

226

=Tg

Appl.

Math.

Modelling,

1996,

Vol. 20, March

of the fluid. If

a2v i au ~1~ a2u arz+vz--+---Q=O

(9)

r

with associated

boundary

conditions

v=o at r=a, r= b, Oszsh and at z = 0, asrsb

(10)

v=rR

on z=h,

boundary and fixed is simpliReynolds to use a

(8)

where p is the constant dynamic viscosity Re = 0 the governing equations reduce to

=-

aU

are

sys-

a+a
I

Additional boundary conditions are required in the gaps _ . but it is noted that if the gaps are zero then an anatytrc solution to the above problem can be obtained (see Appendix). However the solution is very complicated algebraically and it has not been found possible to extend it to allow for nonzero gaps. It is not used in the present work although Turner6 has shown that it gives results in good agreement with those obtained numerically. In this section, standard second-order accurate finite difference representations are used to discretize equation (9), and the resulting

Annular equations are solved iteratively by successive over-relaxation. The solution is obtained on a rectangular mesh with uniform spacing h/N in the z direction and uniform spacing (b - a)/N in the r direction, where N is the number of grid cells in each direction. When the gaps (Y and p are not zero there is a change of type of boundary condition, or possibly a change in gradient of velocity, at the edges of the lid and this leads to singularities in the solution. The treatment of such singularities is discussed by Gladwell and Wait’ and, for the type of singularity encountered in this work, by Fox and Sanka? and Griffiths.’ It is known that the rate of convergence of numerical solutions, as the mesh size is decreased, is significantly reduced when singularities are present, and Griffiths’ for example has shown that it is beneficial to match numerical solutions with analytical solutions valid near the singularity. Four methods are considered for the treatment of the solution near the gaps: (i) ignore the presence of the singularity and impose the boundary condition d~l/dz = 0 in the gaps; (ii) assume that ~1decays linearly to zero across the gaps; (iii) use an approximate analytical solution for u in the gaps; and (iv) match an approximate analytical solution to the finite difference solution near the gaps. Before comparing the solutions obtained by the above methods an approximate solution valid in the gaps is obtained by the method described by Kim.“’ 3.1 Approximate

solution near the edge of the lid

Near the edge of the lid at r = b - p the following variables are introduced: r-b

Xc-

Y=_

z-h

P

Equation

scaled

(11)

P

(9) then becomes

a2 I’

dL>

P

_+_--

p%

b+pXH

ax2

(b+pX)’

+

f37,1 ~

=

8Y2

0

(12) where LJ = c$X, Y). The boundary conditions come, dropping the condition at z = 0, L’ =

(10) be-

\

0

(b-a)

--IX<

at Y=O,

boundary

Y) =

condi-

(16)

where X + iY = cosh( 4 + icp)

(17)

and

From the above solution, expressions for the velocity distribution in the gap and the shear rate on the lid can be found: Ll( r, h) =

z(r,

h) =2R(b-p)(r/h-r)l-pl)

(20) where COS X=(b-r)/P,

(21)

The above equation for av(r, h)/Jz shows clearly the singular nature of the shear rate at the edge of the lid at r = b - p. The torque on the lid is dependent on the integral of the shear rate and is finite since the singularity where d is distance from the edge of is of the form d-‘/’ the lid. It is interesting to note the form that the shear rate takes for very small gap widths. If /3 + 0 then JLl(r, h)/az - l/(b - r), the singularity is of the form d- ‘, and the torque infinite. For nonzero finite (and infiis well nite) gap widths the d- ‘/* type of singularity known for this type of flow.“,‘* Expressions similar to those in equations (19)-(21) can be derived when finding an approximate solution near the inside edge of the lid at r = a + (Y. In method (iii) referred to above, the velocity profiles in the gaps are assumed to be given by equations (19) and (21) and their counterparts for the inside gap. In method (iv), A the expression given by equation (19) for the velocity profile is replaced by (22)

shear rate is then

(23) -1

P For small gap width p, the first approximation (12) is

/ to equation

a2u

(14)

0

-+ay2= 8X2

c=O(b-0)

(14) satisfying

(13)

u=O(b+pX)

u=Oat

u(X,

and the corresponding

-P5yso

and the boundary

The solution to equation tions (15) is

Ll(r, h) =n(b-P)-CpX,

h

at X=0,

CPU

flume flow: P. W. James et al.

conditions

X=0,

(13) are approximated

-m
1

by

(15)

where CD (and C,, its counterpart for the inside gap) is found by matching u(r, z) at a distance of two mesh lengths in the .a direction from the edge of the lid with the numerical solution at that point. The velocity at one mesh length in the z direction below the lid and the velocities in the gap are then found from the adjusted analytic expressions. This is essentially the approach used by Griffiths’ and is illustrated in Figure 2. The results from this method are shown as method (iv)A later. The “exact” values of C, and C, are 4.518 X lop2 and 9.147 X 10m2, respec-

Appl. Math.

Modelling,

1996, Vol. 20, March

227

Annular

flume flow: P. W. James et al. ROTATING

LID

x x I x--SOLUTIONS

Figure 2.

Matching

h

MATCHED

analytic and numerical

HERE

solutions.

values (for N = 40) are tively, whereas the “matched” 4.183 X lo-’ and 9.408 X 10p2, respectively, in ms-‘. 3.2 Comparison

a singularity in Ju(r, h)/az at the edges of the lid of the form log(d) (d being the distance from the edge of the lid), which is weaker than that produced by methods (i), (ii), and (iv). The effect of this singularity on the convergence of the solution would therefore be expected to be weaker too. The integral for r in equation (24) has an infinite integrand at r = a + (Y, r = b - /I and the midpoint rule was used to evaluate the integral rather than, for example, the Trapezium rule or Simpson’s rule. The results obtained for methods (i), (ii), and (iv) A were evaluated this way but use can be made of the known form of the singularity to avoid infinite integrands by writing r in the form

of methods (i)-(iv)

To compare the methods (i)-(iv), attention is focused on evaluation of the torque on the rotating lid, which is given by

r=

ib-P2rrr2rz,(r,

h)dr

(24)

Jata

For simplicity, the gap widths (Y and p are each set to (b - a)/20. For th e miniflumes to be discussed in Section 4, these values are reasonably accurate. The other dimensions required are b = 6.4 cm, a = 2.9 cm, h = 1.6 cm, and the following physical properties are used: p = 884 kg/m3, 0 = 2.308 rad/ set, p = 1.0 Pa . s. These figures are representative of some of the experimental work to be described in Section 4 but are of no special relevance. The geometrical parameters do of course have to be taken into account. In Table 1 we compare the computed values of r for each method for a range of values of N. If it is assumed that the approximate computed torque with a given mesh size, r,, is given by

+ 2?=r( a + cx)‘C,

2r[ r2Tp(r)+ (b-p I

Table 1.

Comparison

where L = (b - p + a + a)/2

20 40 80 160 N+@= 2p

228

Appl. Math. Modelling,

dr

and

= ~Ze(r> h) {T (-

(27)

Tp(r)

= TZs(r7 h)

(28)

The integrands in the above integrals are finite at the edges of the lid so the Trapezium rule can be used. The values of r obtained in this way are presented in Table 1 as (iv)B and are seen to be almost independent of N to an accuracy of two decimal places. No systematic variation with N was found and so extrapolation as N -+ 00 was not attempted. Nevertheless, the values for r obtained by methods (i), (ii), (iv)A, and (iv)B are all within 2% of 8.08 FNm and it is assumed that this is a sufficiently accurate value with which to compare results obtained by HarwellFLOW3D.

Method (i)

(ii)

(iii)

6.090 6.614 7.001 7.289 8.13 1.344

6.090 6.542 6.794 6.927 7.08 1.895

6.090 6.642 7.045 7.335 8.08 1.390

1996, Vol. 20, March

P)2CpdP]

T,(r)

of values of r (in FNm) and 2P for methods (i)-(b).

n

(b -

(26)

(25) where k and p are constant then we can extrapolate and estimate the value of r. as N + m. The values of 2p, based on N = 40, 80, and 160 and of r thus obtained, are also shown in Table 1. If there were no singularities then it would be anticipated that p would be approximately 1. This is clearly not the case and illustrates the effect of the singularities on the rate of convergence of the solutions with decreasing mesh size. In case (ii), p is close to 1 and this is because the change in the boundary conditions at the lid edges involves only a change in the gradient of u(r, h). The changes in the gradients are finite and, using the general solutions obtained by Fox and Sankar,8 lead to

p cash-’

WA 6.452 6.946 7.282 7.507 7.96 1.493

(iv)6 8.066 8.071 8.07 1 8.062 -

Annular Nm (x10-3) 12,

flume flow: P. W. James et al.

N/m-2 i4/

1.7

1.8

2

1.9

2.1

2.2

2.3

2.4

2.5

Volts

Figure 3.

Lid torque - silicone oil.

Figure 5.

The value of r obtained using method (ii) differs significantly from the rest, and it is concluded that the assumption of a linear profile in the gaps is less physically realistic than the other methods. It is interesting to note that method (iv)B gives a result in excellent agreement with the value 8.08 for r using a fairly coarse mesh (N = 20). It is concluded that if the form of the singularity is known then the most efficient method is to match the solution near the singular point(s) to the finite difference solution and to make use of the singular solution when integrating to find the torque. This approach could be used with Harwell-FLOW3D, but if nonzero Reynolds numbers and non-Newtonian fluids are to be considered then it is preferable to devise a grid that is sufficiently refined close to the singularities so that the rapid variation of L’(I, z) near the edges of the lid can be resolved. A grid with cells that contract as we move from the center of the annulus cross-section to the walls has therefore been used. The contraction ratio in the z direction is 1.1 and that in the radial direction is chosen so that the outer edge of the lid coincides with a grid line. With 40 cells in the z direction and 50 in the radial direction a value of r of 7.95 is obtained with Harwell-FLOW3D for the problem considered above. This is considered to be sufficiently close to the value of 8.08 to justify the use of the refined grid

I2

TORSION

BAR

Calibration

curve - Al fluid.

method when using Harwell-FLOW3D mental results.

4. Experimental

to predict experi-

results

4.1 Constant oiscosity silicone oil In Figure 3 the torque on the lid of the miniflume is calculated and compared with experimental results obtained for a silicone oil. The depth of fluid was 1.5 cm, and it is seen that there is generally good agreement between theory and experiment. The oil has an almost constant viscosity of 0.94 Pa s, although rheometer tests show that it is slightly shear thinning. This may account for some of the (slight) discrepancy between theory and experiment. Fine grid solutions of the zero Reynolds number equation and singular solutions of the type presented in Section 3 fall between the Harwell-FLOW3D predictions and the experimental data. 4.2 Shear-thinning

jluids

In previous work3 stress measurements for a shear thinning fluid (polyisobutylene in Dekalin) were made at one point only on the bed of the flume and these were shown to compare well with predicted values. To confirm that the stress distribution over the bed is well predicted, a new annular flume has been constructed in which the “bed” is at the top, the rotating “lid” is at the bottom and bed, and

N/m”2 121

i

101

+

Expe,riment

20

30

---

I

I 40

rm

Figure 4.

Modified

miniature

annular flume.

Figure 6.

Appl.

Bed shear stress at r = 51 mm - Al fluid.

Math.

Modelling,

1996, Vol. 20, March

229

Annular

flume flow: P. W. James et al. Nm (x10-4)

40,

P 15 m 10 I

lf/

5

Nm (xlO”3) 12,

/

E

ri

Ii

OJ

J 10

0

30

20

40

10

0

rm Figure 7.

Figure 9.

the lid and side walls, are free to move relative to one another (Figure 4). The inner and outer radii of the flume and the gaps between the rotating lid and side walls were approximately the same as those given in Figure I, and the experiments were carried out with a 2 cm depth of fluid. Fluid is sealed inside the flume by surrounding the entire apparatus with a reservoir of fluid which is open to the atmosphere. The flume is mounted on a Weissenberg rheogoniometer and a torsion bar is used to measure the torque on the bed. Hot-film shear stress probes are also mounted on the bed to provide simultaneous measurements of stress. The apparatus was used with a shear thinning fluid, a test fluid known as the Al fluid,‘” whose viscosity function was well-fitted by a Cross model:

77”

l+(hj)”

where the constants Q, A, and n are, respectively, 30.45 Pa. s, 4.66 see-‘, and 0.685. The shear stress probes were calibrated in a cone and plate apparatus using a Carri-med rheogoniometer, and the expected linear calibration curve was obtained when y’ was plotted against (voltage)*. It is interesting to note that in this experiment, a plot of shear stress against voltage was almost linear (Figure 5). The results obtained for shear stress, at a point 51 mm from the centre line, and the bed torque are shown in Figures 6 and 7. Excellent agreement between theory and experiment is obtained, confirming that Harwell-FLOW3D is capable of

N/m-2 14

OK 0

2

4

6

8

10

u-m Figure 8.

230

Bed shear stress at r = 51 mm - Sl fluid.

Appl.

Math.

Modelling,

30

40

wm

Bed torque - Al fluid.

T(Y) =

20

1996, Vol. 20, March

Bed torque - Sl fluid.

predicting the bed stress distribution accurately. Note that the lid edge boundary conditions do not significantly affect the predictions of bed stress distribution. For example, in the calculations presented in Section 3 the variation in peak bed stress for methods (i)-(iv) was less than 0.2%. This means that no special treatment or especially refined grid is required to predict the bed stress distribution. Similar calculations were made for another test fluid, known as the Sl fluid, which is again well-fitted by a Cross model but with the parameters r],, = 19.66 Pa s, A = 0.492 set- ’ , n = 0.7. The results are shown in Figures 8 and 9 and it is again seen that excellent agreement is obtained. (Note that in Figure 8 bed stresses were measured over a restricted range of lid speeds because the method of calibration of the shear-stress probes is limited to shear stresses below - 15 Nm-2.) Both of the above fluids are known to be elastic as well as shear-thinning. The Al fluid is highly elastic but the Sl fluid is less so because it is obtained by diluting the Al fluid with a low molecular weight polybutene oil. However as is seen from Figures 6 to 9, ignoring elasticity does not prevent predictions of the bed shear stress from being made that agree with experimental data.

5. Conclusions The main conclusion to be drawn from the above results is that Harwell-FLOW3D can provide good predictions of the bed and lid stress distributions for the slow flow of constant viscosity and shear thinning fluids in annular geometries provided a suitably refined grid is used. Good agreement can also be obtained with a simplified model provided the singularities at the lid edges are taken into account. These results are not immediately applicable to the erosion of cohesive sediments when the overlying fluid is water. In this case the rotational Reynolds numbers are large and the secondary flows are much more significant. However when a cohesive sediment bed erodes, a layer of fluid mud appears immediately above the bed and can be modelled as a shear thinning non-Newtonian fluid. The modification to Harwell-FLOW3D described earlier will then be of direct use.

Annular Acknowledgements The authors have benefitted from the advice of Dr. David Wilton and Dr. David Graham, and D.M.S. is grateful to the University of Plymouth for the award of a research studentship.

flume flow: P. W. James et al.

Turner6 finds approximate values for A, and u(r, z) is evaluated for a number of cases related to the flow in a miniature annular flume. Although the results are generally in excellent agreement with those obtained numerically, the series in (Al) is sometimes very slow to converge.

Appendix When the gap widths CYand p are zero, a solution to the zero rotational Reynolds number problem can be obtained by separation of variables:

cx

u( r, z) =

P,, sinh(A,z)H,(4,r)

n=l

where where

(AlI

sinh( h,h)

A,, is the nth zero of the equation H,( Ar) =.I,( Ar)Y,( Ab) - Y,( Ar)J,(

H,(ha) Ab)

= 0, (A2)

Here .I,,
‘?H,(

A,r)dr (‘43)

”

= /“:{H,( u

A,r)j2dr

and can be reduced to

P, = fi -

“(fnb) [b2J,( A,,b) ?I

- a2J2( A,&]

J’(fnb) [b2Y2( A,b) n

2[J:(A,a)

x -~2YZ(+)l} /i

rr’A;J;(

-J:(A,b)] A\,a)

(A41

References 1. Partheniades, E., Kennedy, J. F., Etter, R. J. and Hoyer, R. P. Investigations of the depositional behaviour of fine cohesive sediments in an annular rotating channel. Hydrodynamics Lab. Rept. 96, Massachusetts Inst. of Tech., Cambridge, MA, 1966 2. Amos, C., Grant, J., Darborn, G. R. and Black, K. Sea carousel-A benthic, annular flume. Es&urine Coastal Shelf&i. 1992, 34, 557577 3. Graham, D. I., James, P. W., Jones, T. E. R., Davies, J. M. and Delo, E. A. Measurement and prediction of surface shear stress in annular flume. AXE J. Hydraul. Res. 1992, 118, 1270-1286 4. Burt, T. N. and Game, A. C. The carousel: commissioning of a circular flume for sediment transport research. Rept. SR 335, Hydraulics Research Ltd., Wallingford, UK, 1985 5. James, P. W. and Jones, T. E. R. The measurement and prediction of surface shear stress induced by Newtonian and non-Newtonian fluids. Proceedings of XIth International Congress on Rheology, Brussels, Elsevier Science Publishers B. V., 1992, pp. 926-928 6. Turner, N. Analytical solution of the problem of laminar flow in an annular flume. Final year honours project report, School of Mathematics and Statistics, University of Plymouth, Plymouth, UK, 1991 7. Gladwell, I. and Wait, R. A Survey of Numerical Methods for Partial Difierential Equations. Clarendon Press, Oxford, 1979, pp. 42-66 in linear elliptic 8. Fox, L. and Sankar, R. Boundary singularities differential equations. J. Inst. Math. Appl. 1969, 5, 340-350 9. Griffiths, D. F. A numerical study of a singular elliptic boundary value problem. J. Inst. Math. Appl. 1977, 19 59-69 10. Kim, M. U. Slow rotation of a disk in a fluid-filled circular cylinder. J. Phys. Sot. Jpn. 1981, 50 4063-4067 11 Levine, H. Skin friction on a strip of finite width moving parallel to its length. J. Fluid Mech. 1957, 3, 145-158 12. Pao, H.-P. Numerical solution of the Navier-Stokes equations for flows in the disk cylinder system. Phys. Fluids 1972, 15, 4-11 13. Hudson, N. E. and Jones, T. E. R. The Al project - An overview. J. Non-Newtonian Fluid Mech. 1993, 46, 69-88

Appl.

Math.

Modelling,

1996, Vol. 20, March

231