Terminal Velocity of Falling Particles

> restart;

The force balance on the particle is

> forcebal := v[t]=sqrt(4*g*(rho[p]-rho)*D[p]/(3*C[D]*rho)):
forcebal;

The Reynolds number is defined for this system by

> reeqn := re = D[p]*v[t]*rho/mu: reeqn;

The parameters in these equations have the following values (in consistent units):

> params := {g=9.81,rho[p]=1800, D[p]=0.208e-3,rho=994.6,mu=8.931e-4};

We substitute these parameters into the force balance to get

> subs(params,forcebal);

We cannot solve this equation directly since the drag coefficient is a function of the terminal velocity through the Reynolds number. We continue by writing a Maple procedure to compute the drag coefficient as a function of the Reynolds number. We have to take this approach since the drag coefficient is computed from different expressions depending on the value of the Reynolds number.

> f := proc(re)
local CD;
if not type(re,numeric) then RETURN(C[D](re)) fi;
if re < 0.1 then
CD := 24/re
elif re >= 0.1 and re <= 1000 then
CD := 24/re*(1+0.14*re^0.7)
elif re > 1000 and re <= 350000 then
CD := 0.44
elif re > 350000 then
CD := 0.19-8e4/re
fi;
end:

Note that we use re for the Reynolds number since Re is a reserved phrase in Maple for the real part of a complex variable. Let us test this proc to see what it does

> f(10);

> f(100);

> f(2000);

> f(re);

The next step is to create a Maple procedure to evaluate the force balance for different values of the terminal velocity. There are many ways to do this; our example follows:

> eqn := proc(vt)
local C, re;
global params;
re:= subs(params,D[p]*vt*rho/mu);
C[D]:=f(re);
subs(params,vt-sqrt(4*g*(rho[p]-rho)*D[p]/(3*C[D]*rho)));
end:

Finally, we invoke Maple's built in floating point solver

> vt1:=fsolve('eqn'(vt),vt=0.00001..0.05);

This is the terminal velocity of the particle under the given conditions. Note that eqn in the call to fsolve is contained within quote marks. If this were not done then Maple would attempt to call the eqn procedure and pass that result to fsolve. Needless to say, that would be a disaster in this case. the quote marks force Maple to re-evaluate the equation every time it tries a new value of the terminal velocity. The units here are m/s. The Reynolds number at this velocity is

> subs(v[t]=vt1,params,reeqn);

and the drag coefficient is

> f(rhs(%));

For part (b) we are asked for the terminal velocity in a centifugal separator where the acceleration is . The only changes needed are to revise the set of parameters (in this case only is changed).

> params := {g=9.81*30,rho[p]=1800, D[p]=0.208e-3,rho=994.6,mu=8.931e-4};

> eqn := proc(vt)
local C, re;
global params;
params := {g=9.81*30,rho[p]=1800, D[p]=0.208e-3,rho=994.6,mu=8.931e-4};
re:= subs(params,D[p]*vt*rho/mu);
C[D]:=f(re);
subs(params,vt-sqrt(4*g*(rho[p]-rho)*D[p]/(3*C[D]*rho)));
end:

Again we call fsolve

> vt2:=fsolve('eqn'(vt),vt=0.0001..1);

The units here are m/s. The Reynolds number at this velocity is

> subs(v[t]=vt2,params,reeqn);

and the drag coefficient is

> f(rhs(%));

>