RK4FlowRates

RK4FlowRates.m

#Flow rate
clear;
clc;

size = 0.00500;
dt = .005;
Volume1 = 120;
Volume2 = 240;

mass = 300;


steps = 1000;

M1 = [];
M2 = [];
P1 = [];
P2 = [];
T = linspace(0,steps*dt,steps);

M1(1) = mass;
M2(1) = 0;
P1(1) = 80000;
P2(1) = 0;
RT = P1(1)*Volume1/(mass);

integrationTemp = zeros(1,4);
flowrate = size*(P1(1)-P2(1))*dt;


method = 1;


for ii = 1:steps
  P1(ii) = M1(ii)*RT/Volume1;
  P2(ii) = M2(ii)*RT/Volume2;
  dP = P1(ii)-P2(ii);

  if method == 1
  #Euler
  flowrate = size*dP*dt;

  elseif method == 2
  #Trapazoid
  integrationTemp(1) = size*dP*dt;
  integrationTemp(2) = size*(dP-integrationTemp(1)/size)*dt;
  flowrate = (integrationTemp(1)+integrationTemp(2))/2;

  elseif method == 3
  #RK4
  integrationTemp(1) = size*dP*dt;
  integrationTemp(2) = integrationTemp(1)*exp(-dt/2);
  integrationTemp(3) = integrationTemp(2)*exp(-dt/2);
  integrationTemp(4) = integrationTemp(3)*exp(-dt);

  flowrate = (integrationTemp(1)+2*integrationTemp(2)+2*integrationTemp(3)+integrationTemp(4))/6;

  elseif method == 4
  flowrate = size*dP*(sinh(dt)-cosh(dt)+1);

  end

  M1(ii+1) = M1(ii) - flowrate;
  M2(ii+1) = M2(ii) + flowrate;
endfor

##figure(1)
##plot(T,M1(1:end-1),T,M2(1:end-1));
##grid on;
##grid minor;
##axis([0,5,0,mass])

figure(2);
plot(T,P1,"g-"); # * o +
grid on;
grid minor;
axis([0,3,0,80000])
mdot = [];
hold on;
title("Analytical");
ylabel("Pressure [Pa]");
xlabel("Time [sec]");

for ii = 1:steps
  mdot(ii) = ((M1(ii)-M1(ii+1)))/dt;
end
##figure(3)
##plot(T,mdot);
##grid on;
##grid minor;