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

In[36]:=

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

Out[36]=

           -7
  1.7017 10   gm

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

In[37]:=

  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}]]

Out[37]=

  5.82866 cm
  ----------
     sec

This would give a Reynolds number of

In[38]:=

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

Out[38]=

  1.93002

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

In[39]:=

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

Out[39]=

  0.234364 cm
  -----------
      sec

and the associated Reynolds number is

In[40]:=

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

Out[40]=

  0.0776038

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

In[41]:=

  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 }

Out[41]=

           0.000553164 cm       0.000553164 cm
  200 cm + -------------- - ----------------------- - 
                  2          (1331.03 Pi t)/sec   2
                Pi          E                   Pi
   
    0.736275 cm t
    -------------
       Pi sec

We wish to find t when y = 0.

In[42]:=

  Solve[height == 0, t]

Out[42]=

                 0.000553164 cm       0.000553164 cm
  Solve[200 cm + -------------- - ----------------------- - 
                        2          (1331.03 Pi t)/sec   2
                      Pi          E                   Pi
   
      0.736275 cm t
      ------------- == 0, t]
         Pi sec

  Solve::tdep: 

     The equations appear to involve transcendental functions of

       the variables in an essentially non-algebraic way.

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.

In[43]:=

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

Out[43]=

        0.000553164      0.000553164      0.736275 t
  200 + ----------- - ----------------- - ----------
              2        1331.03 Pi t   2       Pi
            Pi        E             Pi

In[44]:=

  
  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.

In[45]:=

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

In[46]:=

  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.

In[47]:=

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

Out[47]=

  0.0000105232

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.

In[48]:=

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

Out[48]=

        0.000553164   0.736275 t
  200 + ----------- - ----------
              2           Pi
            Pi

In[49]:=

  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.

In[50]:=

  tsettle = Solve[heightlin == 0, t]

Out[50]=

  {{t -> 853.375}}

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.

In[51]:=

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

Out[51]=

                         -4921
  1.094671782115361521 10

This is certainly a small enough number!!

Up to 6. Applications