#TJHSST ~ Stoke's Law of Settling

stokes.c

#include <stdio.h>
#include <math.h>
#include <gl.h>
#include <device.h>

#define GRAVITY (double)(978.032)   /* centimeters/second^2       */
#define PI (double)(3.141592653589) /* good enough approximation  */
#define DELTAT (double)(0.00025)    /* a very small delta time    */
#define PRECISION (double)(0.00001) /* some level of precision    */
/*  a macro used to calculate the volume of a sphere given its diameter.  */
#define volume(x) ((1.0/6.0)*PI*pow((x),3.0))
#define n20 (double)(1.0019415531)  /* use in tviscosity function */

struct liquid_attr  {
  double viscosity;  /* very viscous fluid is honey, very low is water */
  double density;    /* the mass per unit volume                       */
};

struct sphere_attr  {
  double diameter;   /* it's size                  */
  double volume;     /* volume in cm^3             */
  double mass;       /* amount of material         */
  double density;    /* mass per unit volume       */
  double velocity;   /* velocity of the sphere     */
};

/*  Returns the drag force as calculated by
      Fr = 3*pi*n*v*D
*/
double dragforce(sphere,liquid)
struct sphere_attr sphere;
struct liquid_attr liquid;
{
  return
    3*PI*liquid.viscosity*sphere.velocity*sphere.diameter;
}

/*  Returns the force due to gravity given by
	   1
      Fg = -*pi*(ps - pf)*g*D^3
	   6
*/
double gravforce(sphere,liquid)
struct sphere_attr sphere;
struct liquid_attr liquid;
{
  return (sphere.density-liquid.density)*GRAVITY*sphere.volume;
}

/*  Returns the velocity of the sphere, given by the following
	       F = ma                Newton's Law
     (Fr + Fg)/m = a                 Sum of the forces
	      v1 = v0 + a*DELTAT
*/
double newvelocity(sphere,liquid)
struct sphere_attr sphere;
struct liquid_attr liquid;
{
  double force;

  force=(gravforce(sphere,liquid)-dragforce(sphere,liquid));
  force/=sphere.mass;
  return sphere.velocity+force*DELTAT;
}

/*
  Return the viscosity of water for a given temperature using

  For 0 <= T <= 20
                           1301
     log nT = ---------------------------------------- - 1.30233
              998.333 + 8.1855(T-20.0) + 0.00585(T-20)^2

  For 20 <= T <= 100

         nT   1.3272(20-T) - 0.001053(T-20)^2
    log --- = -------------------------------
        n20              T + 105

        n20 = 1.0019415531
*/
double tviscosity(temperature)
double temperature;
{
  double nT;

  if(temperature<=20.0)  {
    nT=1301.0/(998.333+8.1855*(temperature-20.0)+0.00585*pow(temperature -20.0,2.0));
    nT=nT-1.30233;
    nT=pow(10.0,nT);
  }
  else  {
    nT=(1.3272*(20.0-temperature)-0.001053*pow(temperature-20.0,2))/(temperature+105.0);
    nT=pow(10.0,nT);
    nT=nT*n20;
  }
  return nT;
}

/*
  Return the density of water for a given temperature using

    pT = (999.83952 + 16.945176T - 7.9870401E-3T^2 - 46.170461E-6t^3 +
          105.56302E-9T^4 - 280.54253E-12T^5)/(1 + 16.879850E-3T)
*/
double tdensity(temperature)
double temperature;
{
  double pT;

  pT=(999.83952+16.945176*temperature-
      7.9870401E-3*pow(temperature,2.0)-
      46.170461E-6*pow(temperature,3.0)+
      105.56302E-9*pow(temperature,4.0)-
      280.54253E-12*pow(temperature,5.0))/
      (1+16.879850E-3*temperature);
  return pT/1000.0;
}

/*  Plots the velocity of a particle having the following 
    characteristics.
      ps = density of sphere
       D = diameter of sphere
    Variables based on given parameters are
      pf = density of fluid, based on temperature
       n = viscosity of fluid, based on temperature
       v = velocity, based on previous velocity, and all other
           parameters in previous iterations
    Parameters 
       T = temperature
*/
double tnewvelocity(sphere,liquid,temperature)
struct sphere_attr sphere;
struct liquid_attr liquid;
double temperature;
{
  double force;

  liquid.density=tdensity(temperature);
  liquid.viscosity=tviscosity(temperature);
  force=(gravforce(sphere,liquid)-dragforce(sphere,liquid));
  force/=sphere.mass;
  return sphere.velocity+force*DELTAT;
}

/*  Returns the Reynolds number associated with the given parameters
	   pf*v*D
      Re = ------
	    n 
*/
double reynolds(sphere,liquid)
struct sphere_attr sphere;
struct liquid_attr liquid;
{
  double number;

  number=liquid.density*sphere.velocity*sphere.diameter;
  number/=liquid.viscosity;
  return number;
}

/*  Returns the velocity of a sphere having given 
    characteristics after some period of time.
*/
double plotcourse(density,diameter,temperature)
double density,diameter,temperature;
{
  struct sphere_attr sphere;
  struct liquid_attr liquid;
  double time,velocity;

  sphere.density=density;
  sphere.diameter=diameter;
  sphere.volume=volume(diameter);
  sphere.mass=density*sphere.volume;
  time=0.0;                              /* set initial conds  */
  sphere.velocity=0.0;
  do  {
    velocity=sphere.velocity;
    sphere.velocity=tnewvelocity(sphere,liquid,temperature);
    time+=DELTAT;
  }  while(fabs(velocity-sphere.velocity)>PRECISION);
  return sphere.velocity;
}

main()
{
  int scrx,scry;
  double temperature,diameter,velocity;
  FILE *out;

/*   ginit(); */
/*   color(BLACK); */
/*   clear(); */
/*   cursoff(); */
  out=fopen("what.dat","w");
  for(scry=0;scry<100;scry++)  {
    temperature=scry*1.0;
    for(scrx=1;scrx<26;scrx++)  {
      diameter=scrx/26.0+0.1;
      velocity=plotcourse(1.5,diameter,temperature);
      fprintf(out,"%d %d %10.8lf\n",scrx,scry,velocity);
/*       mapcolor(10,(int)(temperature*2.0),128,(int)(velocity)); */
    }
  }
/*   curson(); */
  fclose(out);
}

A C program found on an old 5.25 floppy disk from when I was a nerd at Thomas Jefferson High School for Science and Technology between 1988 and 1992 for a lab in my Geoscience class.

The program calculates (and displays) various parameters related to "Stoke's Law of Settling."

http://en.wikipedia.org/wiki/Stokes%27_law

Returns the drag force as calculated by

Fr = 3*pi*n*v*D

http://en.wikipedia.org/wiki/Stokes%27_law#Terminal_velocity_and_settling_time

Returns the force due to gravity given by

     1
Fg = -*pi*(ps - pf)*g*D^3
     6

http://en.wikipedia.org/wiki/Reynolds_number

Returns the Reynolds number associated with the given parameters

     pf*v*D
Re = ------
       n

Return the viscosity of water for a given temperature using

For 0 <= T <= 20

                       1301
log nT = ---------------------------------------- - 1.30233
          998.333 + 8.1855(T-20.0) + 0.00585(T-20)^2

For 20 <= T <= 100

     nT   1.3272(20-T) - 0.001053(T-20)^2
log --- = -------------------------------
    n20              T + 105

    n20 = 1.0019415531

Return the density of water for a given temperature using

pT = (999.83952 + 16.945176T - 7.9870401E-3T^2 - 46.170461E-6t^3 +
      105.56302E-9T^4 - 280.54253E-12T^5)/(1 + 16.879850E-3T)

Anyone can do whatever they'd like to with this program--if anything.

I remain a science nerd.