Solved: MATLAB While Loop Help for Final Time, Mass & Temp

[Juanka]

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

 

[Mark44]

I'm not a Matlab expert, but here are a couple of things I noticed. I also fixed up your indentation to make the code easier to read.

while P2<P1-.01
    t==0+dt; % [color="red"]Should be t = 0 + dt; or better, t = dt;[/color]
    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; % [color="red"]Remove the ;   Also, could change to if Prat < Pcrit[/color]
            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

 

[FourierFaux]

I suspect the issue is in your first while loop.

Your t array keeps getting redefined from:

t(n) = n*dt; % in your second loop.

to

t == 0 + dt; % '==' means defined to be. This will spit out a boolean: True or False

You're graphing all of the variables with respect to t. (I suspect that you're destroying your t array every time you exit out of the second loop and go back into your first)

And... it looks like you found an infinite series solution to a differential equation. (Did you do all of that paperwork by hand?)

I think everything else should work properly... (I haven't tested it personally, so I may be wrong). You can test to see if everything is working properly by inserting little variable values in the code and see if they behave the way you expect them to.

Out of curiosity, (this has nothing to do with your code) can I see the differential equation you're trying to solve?

 

[Juanka]

I have gone away with my first time step, I am just leaving the second I got confused when coding and didn't even realize I had coded it twice.

I do have the work done on paper by hand, and I had already done what you suggested for checking my program and the first iteration that my program runs gives me the same values I had on paper.

As far as solving a differential equation, I am using the flow equation below (which in the code I sent you had an error [ i had a times *18 when it should have been 144 (the conversion from psi to psf ]) as well as the energy equation used in the program to find the time, mass, and final temperature of a pipe filling a tank reservoir of a volume of 70ft^3.

However my answers are still wrong, I am more familiar with the thermodynamics behind this problem than the coding. I really do not know where the problem is. I am going to attach the problem statement so that someone may better understand what is going on and possibly help me find the error.

Problem Statement

A 70 ft3 rigid insulated tank contains air at 14.7 psia and 80 º F. The tank is connected to a supply line through a valve. Air is flowing in the supply line at 75 psia and 80 º F. The valve is opened, and air is allowed to enter the tank until the pressure in the tank reaches the line pressure, at which point the valve is closed. The diameter of the connecting pipe is 1 inch.

Write a computer program to model the pressurization process in the tank. The computer program will calculate the pressure, temperature, mass flowrate into the tank and resident mass during the process. Determine the final temperature, mass of air in the tank and time required to pressurize the tank.


However I noticed I defined Prat and Pcrit outside of my loop is this restricting the values of my flow equation just to be the initial conditions or should it still be changing with each iteration?

I am re posted my updated code so everything can be seen visually

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

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

    error=1;

       n=n+1;

    while error>0.001

        %Equations

       mdot=A1*sqrt(((2*k)/(k-1))*P1*144*g*rho1*(Prat)^(2/k)*(1-(Prat))^(k-1/k));%mass flow

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

 

[Juanka]

I have gone away with my first time step, I am just leaving the second I got confused when coding and didn't even realize I had coded it twice.


However my answers are still wrong, I am more familiar with the thermodynamics behind this problem than the coding. I really do not know where the problem is. I am going to attach the problem statement so that someone may better understand what is going on and possibly help me find the error.

Problem Statement

A 70 ft3 rigid insulated tank contains air at 14.7 psia and 80 º F. The tank is connected to a supply line through a valve. Air is flowing in the supply line at 75 psia and 80 º F. The valve is opened, and air is allowed to enter the tank until the pressure in the tank reaches the line pressure, at which point the valve is closed. The diameter of the connecting pipe is 1 inch.

Write a computer program to model the pressurization process in the tank. The computer program will calculate the pressure, temperature, mass flowrate into the tank and resident mass during the process. Determine the final temperature, mass of air in the tank and time required to pressurize the tank.



However I noticed I defined Prat and Pcrit outside of my loop is this restricting the values of my flow equation just to be the initial conditions or should it still be changing with each iteration?

I am re posted my updated code so everything can be seen visually

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

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

    error=1;

       n=n+1;

    while error>0.001

        %Equations

       mdot=A1*sqrt(((2*k)/(k-1))*P1*144*g*rho1*(Prat)^(2/k)*(1-(Prat))^(k-1/k));%mass flow

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

 

[FourierFaux]

Perhaps if you can walk me through your thinking in your solution to this problem. The issue with trying to program something like this in Matlab, when you haven't really done problems like this before, is that when you tell Matlab to do something, it does EXACTLY what you tell it to.

The problem is: If you're not quite experienced with trying to translate your model into code, you may not be translating it correctly. So, if you wouldn't mind, please walk me through the mathematics behind your model.

Are you using the fluid equations with a temperature diffusion term? Or are you using something simpler?

 

[FourierFaux]

I found the paper you're basing this off of! I'll get back to you. This seems like an interesting little toy problem. :)

 

[FourierFaux]

To answer your question, now that I feel that I understand what you're doing. Yes. You're P1 (which represents the pressure in your chamber) should be evolving in time and changing your flowrate. The way you have it written, the flowrate is going to be constant. Therefore it will take only 5.38 seconds for you're reservoir to fill instead of the 8 seconds that you say it should. (I copied your code and have played with it a bit, I hope that's ok.)

 

[FourierFaux]

This is 5.38 seconds with a constant flowrate.

If the flowrate is slowing down, time should be more than that.

Do you know how to do this in your code?

 

[Juanka]

No, I do not that is where I am stuck at, the P1 should be constant at 75 psi, the P2 is what should change with time, and yes I encourage you to play with my code, however the final time should be around 15 seconds*** with a mass of about 20 lbs and a temperature around 700 Rankine. The mass and temperature I get are very similar but my time is way off (5.34 seconds as you said)

 

[FourierFaux]

Your Mass Flow Rate:
mdot=A1*sqrt(((2*k)/(k-1))*P1*144*g*rho1*(Prat)^(2/k)*(1-(Prat))^(k-1/k));%mass flow


What I suspect it should be:
mdot=A1*sqrt(((2*k)/(k-1))*P1*144*g*rho1*(P2/P1)^(2/k)*(1-(P2/P1))^(k-1/k));%mass flow


Another Observation:
The mass flow rate in equation 2 of the paper doesn't have this Area term that your equation seems to. Why is this?

 

[Juanka]

The mass flow in the paper is Mdot/A, the paper you found online Ibelieve was written by my professor (Majumdar), but in this example we have to multiple the flow by A to find the mass flow. However, I have been trying my code as u just suggested and the answer still is about the same, are you getting the same results?

 

[FourierFaux]

Yep. It means that there is more that we're missing. I'm getting a new time of 6.33 seconds... better, but not the 20 you mentioned.

 

[Juanka]

Yes I am getting 6.33 with the values you said, I believe it is because the equation is not converging at each time interval but I am not familiar enough with MATLAB to modify my program to incorporate this.

 

[Juanka]

It looks like the flow equation is still not converging with time even using the flow equation in the form of

mdot=A1*sqrt(((2*k)/(k-1))*P1*144*g*rho1*(P2/P1)^(2/k)*(1-(P2/P1))^(k-1/k));%mass flow

However if you substitute an alternate form of a flow equation

mdot=Cd*A1*sqrt(2*rho1*g*(P1-P2)*144); %flow equation

The graph displays a plot characteristic to that from the report posted online, however, we have to use the same flow equation as in the report, I do not know what is going wrong here.

 

[FourierFaux]

I've solved DE's numerically on Matlab before; from reading the paper, it's not quite clear to me what the convergence algorithm should look like. He's described the other equations well enough. It would be really useful to find that copy of "Unsteady Flow Thermofluid Mechanics" by Moody that he's basing the formulae on, I'd imagine that the model is explained in more detail there.