# Plotting Terminal velocity Vs particle diameter

4 views (last 30 days)
Ahmed Alshehhi on 3 Dec 2019
Commented: Ahmed Alshehhi on 3 Dec 2019
I am trying to plot Terminal velocity Vs particle diameter put the plot shows nothing
%Calculation of terminal velocity
g= 9.8; % in m/s^2
rau_p= 7850; %particle density in Kg/m^3
rau_f= 1000; % fluid density in Kg/m^3
vis_f= .001; % fluid viscoisty in Pa.s
C_D=zeros(1); %initial guess
for Dp=0.0001:.00001:0.1 %in m
Vt=sqrt(4*g*(rau_p-rau_f)*Dp/(3*C_D*rau_f)); %terminal velocity in m/s
Re=rau_f*Vt*Dp/vis_f;
if Re<=.1
C_D=24/Re;
Vt=sqrt(4*g*(rau_p-rau_f)*Dp/(3*C_D*rau_f))
elseif Re>0.1 && Re<=1000
C_D= (24/Re)*(1+.14*(Re^.7));
Vt=sqrt(4*g*(rau_p-rau_f)*Dp/(3*C_D*rau_f))
elseif Re>1000 && Re<=350000
C_D=0.445;
Vt=sqrt(4*g*(rau_p-rau_f)*Dp/(3*C_D*rau_f))
elseif Re>350000 && Re<=1000000
C_D=interp1([350000 400000 500000 700000 1000000],[.396 .0891 .0799 .0945 .11],Re);
Vt=sqrt(4*g*(rau_p-rau_f)*Dp/(3*C_D*rau_f))
else %Re>1000000
C_D= .19-((8*10^4)/Re);
Vt=sqrt(4*g*(rau_p-rau_f)*Dp/(3*C_D*rau_f))
end
end
x=Dp
y=Vt
plot(x,y)
xlabel('D(m)')
ylabel('Vt(m/s)')

Stijn Haenen on 3 Dec 2019
With this script you overwrite vt and dp every loop, so at the end you only have 1 single value for vt and dp and you cannot plot a single value.
I have added line 6,7 28,29 and 33 to get a working script (see below):
g= 9.8; % in m/s^2
rau_p= 7850; %particle density in Kg/m^3
rau_f= 1000; % fluid density in Kg/m^3
vis_f= .001; % fluid viscoisty in Pa.s
C_D=zeros(1); %initial guess
plotdata_x=[];
plotdata_y=[];
for Dp=0.0001:.00001:0.1 %in m
Vt=sqrt(4*g*(rau_p-rau_f)*Dp/(3*C_D*rau_f)); %terminal velocity in m/s
Re=rau_f*Vt*Dp/vis_f;
if Re<=.1
C_D=24/Re;
Vt=sqrt(4*g*(rau_p-rau_f)*Dp/(3*C_D*rau_f))
elseif Re>0.1 && Re<=1000
C_D= (24/Re)*(1+.14*(Re^.7));
Vt=sqrt(4*g*(rau_p-rau_f)*Dp/(3*C_D*rau_f))
elseif Re>1000 && Re<=350000
C_D=0.445;
Vt=sqrt(4*g*(rau_p-rau_f)*Dp/(3*C_D*rau_f))
elseif Re>350000 && Re<=1000000
C_D=interp1([350000 400000 500000 700000 1000000],[.396 .0891 .0799 .0945 .11],Re);
Vt=sqrt(4*g*(rau_p-rau_f)*Dp/(3*C_D*rau_f))
else %Re>1000000
C_D= .19-((8*10^4)/Re);
Vt=sqrt(4*g*(rau_p-rau_f)*Dp/(3*C_D*rau_f))
end
plotdata_x(end+1)=Dp;
plotdata_y(end+1)=Vt;
end
plot(plotdata_x,plotdata_y)
xlabel('D(m)')
ylabel('Vt(m/s)')

#### 1 Comment

Ahmed Alshehhi on 3 Dec 2019
It worked thanks alot !