On Sun, 4 Mar 2001 02:42:20 +0100, "Rune" <run### [at] inamecom>
wrote:
>line. I'm can't figure out if it's supposed to behave that way or not? How
>should the particles move when wind is used?

The drag coefficient is highly speed dependent. A side wind on a
falling sphere will initially give a large wind force, then decrease
as the sphere accelerates to steady state.

Below I've included a drag coefficient function based on experimental
data. I don't have the time to convert it to pov syntax, but you
probably know what to do. (CD is the output, so when this is written
as a macro just type CD at the end.) Input is the Reynold's number. If
hydrodynamic mass and buoyancy is neglected, the following
differential equation is obtained:
dV/dt=-g*y-C*Vr*abs(Vr)*CD
where 
V is speed <vx,vy,vz>
y=<0,1,0>
C=0.75*rho/d
rho=rho_air/rho_sphere
Vr=V_sphere-V_air
abs(V)=vlength(V)
CD=CDkule(Re)=CDkule(abs(Vr)*d/nu)
d is sphere diameter and nu is kinematic viscosity i.e 1.5e-5 for air
when everything is in si units. 

The differential equation can be solved within pov using Euler, or
even better, Heun's method. If you're not familiar with numerical
integration maybe someone else here can help you. I don't have the
time at the moment.

sig

Short explanation of the syntax below:
array=[43 23 13 4]; // defines an array (or row matrix actually)
array(3) // gives 13


%======================== CDkule ==========================
% Computes the drag-coefficient of a sphere as a function
% of the Reynolds nunber Re. Curve fitted after fig. A-56
% in "Fluid Mechanics & Hydraulics", Schaum's Solved Problems
%=============================================================


%-----------------------
function val = qparab(y,x) 
val = y(1)+x*(y(2)+(y(3)+ y(4)*x)*x);


function CD = CDkule(Re)
a2=[0.658 34.87 -14.05 4.22];
a3=[2.41 -17.08 43.72 -30.41];
a4=[11.932 -8.472 2.031 -0.1584];
a5=[14.93 -7.332 1.1905 -0.06335];
if (Re<=0.0)
   CD = 0.0;
   return
end
if (Re> 8.0e6) 
   CD=0.2;
   return
end
x0 = log(Re)/log(10.0);
if (Re> 0.0 & Re<=0.5)
   CD=24.0/Re;
   return
end
if (Re> 0.5 & Re<=100.0)
    CD=qparab(a2,1.0/Re);
    return
end
if (Re> 100.0 & Re<=1.0e4)
    CD=qparab(a3,1.0/x0);
    return
end
if (Re> 1.0e4 & Re<=3.35e5)
    CD=qparab(a4,x0);
    return
end
if (Re> 3.35e5 & Re <=5.0e5)
    x1 = log(Re/4.5e5)/log(10.0);
    CD=91.08*x1^4 + 0.0764;
    return
end
if (Re> 5.0e5 & Re<=8.0e6)
    CD=qparab(a5,x0);
    return
end


-- 
ICQ 74734588

Post a reply to this message