[cenker2018]

Imagine a steel sphere of diameter d is dropped in air with density Pf = 1.29 kg/ m^3 and kinematic viscosity coefficient v = l .49xl0^-5 m^2/s. The density of the steel sphere is p = 8xl0^3 kg/m^3 and the gravitational acceleration is g = 9·.81 m/s^2 .

Now, write and run a matlab code by employing the fourth order Runge-Kutta method in order to examine numerically the motion of the steel sphere. Run the code until Tmax = 5 seconds is reached, h = 0.1 s. is suggested for the time step. Plot the time history of the displacement and the velocity for steel spheres of different diameters; d = 0.07, 0.02, 0.01 and 0.001 m.

WHAT IS THE WRONG WITH THE CODES ?

% compare the motion of a steel sphere dropped in air with that in vacuum,
% based on fourth-order runge-kutta method

---main function---------
% specify data
t0=0;
z0=0;
v0=0;
rho=8000;
rhof=1.22;
g=9.8;
h=0.1;
tmax=10;
rhobar=rhof/rho;
nu=0.0000149;
a=1+rhobar/2;
b=(1-rhobar)*g;
d=0.01;
c=3*rhobar/(4*d);

% initialize the problem
t=t0;
z=z0;
v=v0;
re=v*d/nu;

% compute position and velocity in vacuum, and print data
n=1;

while t<=tmax
zv=z0+v0*t+g*t*t/2;
vv=v0+g*t;

% runge-kutta method
d1z=h*v;
d1v=h*f(v,a,b,c,d,nu);
d2z=h*(v+d1v/2);
d2v=h*f(v+d1v/2,a,b,c,d,nu);
d3z=h*(v+d2v/2);
d3v=h*f(v+d2v/2,a,b,c,d,nu);
d4z=h*(v+d3v);
d4v=h*f(v+d3v,a,b,c,d,nu);
t=t+h;
n=n+1;

z=z+(d1z+2*d2z+2*d3z+d4z)/6;
v=v+(d1v+2*d2v+2*d3v+d4v)/6;
re=v*d/nu;
end

------function-------
function f=f(v,a,b,c,d,nu)
% f=function on right-hand side of equation (1.2.3)
% w=velocity
% r=reynolds number

r=abs(v)*d/nu;

% stokes formula cannot be used when r=0. in this case the value zero is
% assigned to cd
if r==0
cd=0;
end

if r>0 && r<=1
cd=24/r;
end

if r>1 && r<=400
cd=24/r^0.646;
end

if r>400 && r<=3e5
cd=0.5;
end

if r>3e5 && r<=2e6
cd=3.66e-4*r^0.4275;
end

if r>2e6
cd=0.18;
end

f=(b-c*v*abs(v)*cd)/a;