Time step help within a while loop

Ok I have the following matlab code I have written I am trying to find the time,mass, and final temperature of a tank being filled by a supply line, however I need help defining my time step in matlab. or I need help with the while loop itself the results I am supposed to be getting are Final Time= about 15 seconds Final Temp= about 700 R Final Mass= about 20 lbs %%%%%%%%%%%%%% My program is getting these answers but it is not converging at each time step so my plots are all linear when the should not be. Please any suggestions will be greatly appreciated!!

P1=75; % in psi

T1=(80+459.67); %in R

A1=(pi*1)/(4*144); % Area in ft^2

D1=1; % Diameter in inch

P2=(14.7); % in psi

T2=T1; %in R

V2=70; % in ft^3

Cp=0.24; % in Btu/lb-R

Cv=0.17; % in Btu/lb-R

Cd=0.6; % is unitless

R=53.33; % in ft-lb/lb/R

rho1=(P1*144)/(R*T1); % density in lb/ft^3

rho2=(P2*144)/(R*T2); % density in lb/ft^3

mass2=(rho2*V2); % mass in lb

g=32.174; % in ft/sec^2

dt=0.01; % change in time intervals
t=0;
n=1;
k=1.4;
m(1)=mass2;

P(1)=P2;

T(1)=T2;
Pcrit=(2/(k+1))^(k/(k-1));
Prat= (P2/P1);

while P2<P1-.01
    t==0+dt;

      error=1;
         n=n+1;
      while error>0.001
          %Equations
  mdot=A1*sqrt(((2*k)/(k-1))*P1*18*g*rho1*Prat^(2/k)*(1-(Prat))^(k-1/k));
if Prat <= Pcrit;
            Prat= Pcrit;
        end   
mass2new=mass2+mdot*dt;

          u=(mass2*Cv*T2+mdot*Cp*T1*dt)/mass2new; %energy equation
          T2=u/Cv;
          P2new=mass2new*R*T2/(V2*144);
          error=abs(P2new-P2)/P2;
          P2=P2new;
      end
      mass2=mass2new;
      T(n)=T2;
      P(n)=P2;
      m(n)=mass2new;
      t(n)=n*dt;
end

%%Output

fprintf('The final temperature = %7.3f R\n',T2)

fprintf('The mass of air in tank = %7.3f lb\n',mass2new)

fprintf('The time required to pressurize the tank = %7.3f s\n',t(n))

figure(1)

plot(t,T,'g','Linewidth',2)

grid

ylabel('Temperature (R)')

xlabel('Time (s)')

title('Time Vs Temperature')

figure(2)

plot(t,P,'b','Linewidth',2)

grid

xlabel('Time (s)')

ylabel('Pressure (psi)')

title('Time Vs Pressure')

figure(3)

plot(t,m,'r','Linewidth',2)

grid

xlabel('Time (s)')

ylabel('Mass (lb)')

title('Time Vs Mass')

I am getting the correct answers, but my plots are wrong because my while loop is not converging at each time step. Please help!

Answer 1

[the cyclist]

I wonder if the location of your assignment

 Prat= (P2/P1);

is incorrect. Should that ratio stayed fixed for the entire program?

It was also not clear to me if your P2=Pnew statement is in the right place.

[John]

I think you are correct about the Prat, since I am defining it outside the loop it is not taking in to affect the new pressures correct? and do I just need to redefine the location of Prat=(P2/P1) to be in the while loop?

[the cyclist]

I have to admit that I have not spent the time to understand your problem deeply. I don't really know exactly when Prat should be changing.

[John]

The pressure ratio should be less than the critical pressure. If not the assumption should be made that the pressure ratio is equal to critical pressure. I have critical pressure (Pcrit) already defined in the program.