Matlab- Thermodynamic filling of tank.

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

k=1.4;% gamma

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;

m(1)=mass2;

P(1)=P2;

T(1)=T2;


while P2<P1-.01

    error=1;

       n=n+1;

    while error>0.001

        %Equations

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

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