MATLAB: Plotting Terminal velocity Vs particle diameter

plottingvelocity

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)')

Best Answer

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)')
Related Question