Example: Sand settling in water

Suppose the Toronto water authority uses Lake Ontario as a water source for the city. After collecting the water, they must allow sand particles to settle. Our problem is to find out how long it will take the sand to settle in holding tanks that are 200 cm deep. The diameter of fine-grained sand is about 5 10^-3 cm and its density is 2.6 gm/cm3. We must determine if we can use the linear or quadratic drag model by checking the size of the Reynolds number. Unfortunately, we must first know the velocity, since Re = r v D/m. But we need to know which model to use before we can find the velocity! We will resolve this dilemma by approximating the velocity with the terminal velocity for each model. Then we can check to see if the terminal velocity for either model is consistent with its valid range of Reynolds numbers. We need to calculate the mass of a sand particle in order to find its terminal velocity. Multiplying the density by the volume, we get


  Clear[cm,gm]
  dsand = 5.0 10^-3 cm;
  msand = N[(2.6 gm/cm^3) * (1/6 Pi dsand^3)]

For the quadratic drag model, using CD = 0.5, the terminalvelocity is


  Clear[cm,sec]
  g = 980 cm/sec^2;
  vtsand =N[PowerExpand[Sqrt[m g/k2] 
              /. {m -> msand, rho -> 1.00 gm/cm^3, d -> dsand}]]

This would give a Reynolds number of


  Clear[cm,gm,sec]
  rho vtsand d/mu /.
           {rho -> 1.00 gm/cm^3, d -> dsand,
            mu -> 1.51 10^-2 gm/(cm sec)}

This is below the appropriate range, 10^3 < Re < 10^5, for the quadratic drag model. Since the velocity of the sand particles will always below the terminal velocity, we must reject the quadratic drag model for this problem.

Repeating this procedure for the linear drag model, we get


  vtsand = N[m g/k1 /. {m -> msand, d -> dsand, 
                      mu -> 1.51 10^-2 gm/(cm sec)}]

and the associated Reynolds number is


  rho vtsand d/mu /.
           {rho -> 1.00 gm/cm^3, d -> dsand,
            mu -> 1.51 10^-2 gm/(cm sec)}

This is well within the range of valid Reynolds numbers, Re < 0.5, for the linear drag model.

Our problem asks when a particle reaches the bottom of the tank; thus we want to find the value of t when y = 0. Substituting v0 = 0 cm/sec and y0 = 200 cm into the solution from Section 5.1, we get


  Clear[k,m,d,g,mu,v0,y0,gm,cm,sec]
  height1a = height1 /. k -> k1;
  height = height1a /. { m -> msand, d -> dsand,
              g -> 980 cm/sec^2, mu -> 1.51 10^-2 gm/(cm sec),
              v0 -> 0 cm/sec, y0 -> 200 cm }

We wish to find t when y = 0.


  Solve[height == 0, t]

Mathematica cannot find the exact solution because t appears both linearly and as an exponent in the expression for the height. We must find an approximate solution. Let's plot the graph in order to get an initial estimate.


  
  heightdimless = Simplify[height/cm /. t -> t sec]


  
  heightplot = Plot[heightdimless, {t,0,1000}, 
     PlotRange -> {0,200}];

Hey, it looks like a linear function! What's going on here? Let's examine the formula for the height more carefully. The first two terms are constants and the fourth term is linear in t. The denominator of the third term is an exponential with a large coefficient of t. This exponential grows very quickly and thus the third term rapidly approaches zero. Let's look at the height at small values of t, to see how quickly the exponential approaches zero.


  Plot[heightdimless, {t ,0,1}];


  Plot[heightdimless, {t,0,0.001}, PlotRange -> {199.9999,200}];

It appears that the graph becomes linear at about t = 0.0004. We can verify this numerically.


  expterm = 
     0.0005531636448519221/(E^(1331.02503176545 Pi t) Pi^2);
  N[expterm /. t -> 0.0004]

Try other values of t. You will find that t could be chosen to be even smaller than 0.0004. Physically, this means that the sand particles approach their terminal velocity very quickly.

The exponential term becomes insignificant for such small values of t that we can safely ignore this term. Let's plot our linear approximation and the the actual height on the same graph to verify the approximation.


  heightlin = heightdimless[[1]] + heightdimless[[2]] +
              heightdimless[[4]]


  heightlinplot = Plot[heightlin, {t,0,1000},
       PlotRange -> {0,200}, DisplayFunction -> Identity,
       PlotStyle -> {{RGBColor[0,0,1],Thickness[0.01]}}];
  Show[{heightplot, heightlinplot}, 
       DisplayFunction -> $DisplayFunction];

Yes, we can safely use the linear approximation. Let's try to find the settling time again.


  tsettle = Solve[heightlin == 0, t]

The sand will settle in 853 seconds. we can find the value of the exponential term at t = 853 to verify that it is negligible.


  N[expterm /. t -> tsettle[[1,1,2]]]

This is certainly a small enough number!!

Up to 6. Applications