MATLAB: The problem is the graph comes up as blank

graphplot

g= 9.8; % in m/s^2
rop= 7850 %particle density in Kg/m^3
rof= 1000; % fluid density in Kg/m^3
vis= 0.001; % fluid viscoisty in Pa.s
cd=zeros(1); %initial guess
for Dp=0.0001:0.0002:0.150 %mm
ut=sqrt((4*g*(rop-rof)*Dp)./(3*cd*rof)); %terminal velocity in m/s
Re=rof*ut*Dp/vis;
if Re<=0.1
cd=24./Re;
elseif (Re>0.1) & (Re<=1000)
cd=(24/Re)*(1+0.14*(Re^.7));
elseif (Re>1000) & (Re<=350000)
cd=0.445;
elseif (Re>350000) & (Re<=1000000)
cd=(-2*10^7)*re+0.2871;
else %Re>1000000
cd=0.19-((8*10^4)/Re);
end
end
figure
plot(Dp,ut)
xlabel('Dp(m)')
ylabel('ut(m/s)')

Best Answer

Initialise ‘cd’ with eps instead of 0, subscript ‘cd’ and ‘ut’, iterate over ‘k’ and not ‘Dp’ (since you need to store the values of the variables), and it works: 
g= 9.8; % in m/s^2
rop=  7850; %particle density in Kg/m^3
rof= 1000; % fluid density in Kg/m^3
vis= 0.001; % fluid viscoisty in Pa.s
Dpv=0.0001:0.0002:0.150; %mm
cd=zeros(size(Dpv))+eps; %initial guess
for k = 1:numel(Dpv)
    Dp = Dpv(k);
    ut(k)=sqrt((4*g*(rop-rof)*Dp)./(3*cd(k)*rof)); %terminal velocity in m/s
    Re=rof*ut(k)*Dp/vis;
    if Re<=0.1 
        cd(k)=24./Re; 
    elseif (Re>0.1) & (Re<=1000)
        cd(k)=(24/Re)*(1+0.14*(Re^.7));
    elseif (Re>1000) & (Re<=350000)
        cd(k)=0.445;
    elseif (Re>350000) & (Re<=1000000)
        cd(k)=(-2*10^7)*Re+0.2871;
    else %Re>1000000
        cd(k)=0.19-((8*10^4)/Re); 
    end
end
figure
 plot(Dpv,ut)
xlabel('Dp(m)')
ylabel('ut(m/s)')
.
Related Question