[Mateja]

I am trying to simulate Sod's shock tube using PISO and FVM in MATLAB. I tried making this code but the results are wrong. Can someone check the code and let me know where the error is
Thank you.

% Code for Inviscid flow in a shock tube
close all
clear all
clc
n_points =101; % These are the velocity nodes P nodes would be n_points+1
dom_length = 0.35;
dx = dom_length/(n_points-1);
x = 0:dx:dom_length; %These are the velocity nodesor cell faces
x_Press(1:n_points-1) = x(1:n_points-1)+(dx/2);

% initialize the variables
dt = 2.05*10.^-6;
% dt = 1;
time = 0;
final_time = 1;
highP = 6.897*10.^4;
lowP = 6.897*10.^3;
highTemp = 288.89;
lowTemp = 231.11;
% high_Density=1
% low_Density = 0.125
Gas_Constant = 287.049;
gamma = 1.4;
guess_Vel(1:n_points) = 0;

for index = 1:n_points-1
if x_Press(index) < 1.524*10.^-1
g_Pressure(index) = highP;
else
g_Pressure(index) = lowP;
end
end
Ini_P = g_Pressure;
% mean = mean(g_Pressure)

% Specify the initial conditions
for index = 1:n_points-1
if x_Press(index) < 1.524*10.^-1
g_Pressure(index) = g_Pressure(index);
Temperature(index)= highTemp;
g_Density(index) = (g_Pressure(index))./(Gas_Constant*Temperature(index));
Sound_Speed(index)= sqrt(gamma*Gas_Constant*Temperature(index));
else
g_Pressure(index) = g_Pressure(index);
Temperature(index)= lowTemp;
g_Density(index) = (g_Pressure(index))./(Gas_Constant*Temperature(index));
Sound_Speed(index)= sqrt(gamma*Gas_Constant*Temperature(index));
end
end

Ini_P = g_Pressure;
while time < final_time
%Dirichlet BC on the 1st and last velocity bound node while the
%momentum equation will be solved for 2nd to 2nd last velocity CV
% Find density on velocity nodes for the transient term
for index = 1:n_points
if index ==1
Den_on_Vnodes(index)= highP/(Gas_Constant*highTemp);
elseif index == n_points
Den_on_Vnodes(index)= lowP/(Gas_Constant*lowTemp);
else
Den_on_Vnodes(index) = (g_Density(index-1)+g_Density(index))/2;
end
end

% Find velocity on scalar nodes for the non linear term
for index = 1:n_points-1
V_On_Scalars(index) = (guess_Vel(index) + guess_Vel(index+1))/2;
Den_Vel(index)= g_Density(index)*V_On_Scalars(index);
end

% Find gradient of pressure on velocity nodes
for index = 1:n_points
if index == 1
gradP(index) = 0;
elseif index ==n_points
gradP(index) = 0;
else
gradP(index) = (g_Pressure(index-1) - g_Pressure(index));
end
end

% Combine RHS of the x mom eq
for index = 1:n_points
RHS(index)= gradP(index) + (Den_on_Vnodes(index)*guess_Vel(index)*(dx/dt));
end

% complete the equations
V_Pred= sym('V',[1 n_points]);
for index = 2: n_points-1
if index ==2
eqn(index) = ((Den_on_Vnodes(index)*(dx/dt))+Den_Vel(index))*V_Pred(index) == RHS(index);
else
eqn(index) = -(Den_Vel(index-1)*V_Pred(index-1))+(((Den_on_Vnodes(index)*(dx/dt))+Den_Vel(index))*V_Pred(index)) == RHS(index);
end
end

[A,B] = equationsToMatrix(eqn(2:n_points-1),V_Pred(2:n_points-1));

Pred_Velo(2:n_points-1) = linsolve(A,B);
Pred_Velo(1) = 0;
Pred_Velo(n_points) = 0;
Pred_Velo = Pred_Velo';

% 1st pressure correction
% find new density and corr density
for index = 1:n_points-1
DenS(index) = g_Pressure(index)/(Gas_Constant*Temperature(index));
CorrDen(index) = ((DenS(index) - g_Density(index))*dx)/dt;
end

for index = 1:n_points-1
P_Source(index) = -CorrDen(index)+ (Den_on_Vnodes(index)*Pred_Velo(index)) - (Den_on_Vnodes(index+1)*Pred_Velo(index+1));
end

% complete the equations
P_Corr= sym('P',[1 n_points-1]);
for index = 1: n_points-1
if index == 1
P_eqn(index) = 2*P_Corr(index) - P_Corr(index+1) == highP+ 2*g_Pressure(index) - g_Pressure(index+1)- highP + P_Source(index)*(dx/dt);
elseif index == n_points-1
P_eqn(index) = 2*P_Corr(index) - P_Corr(index-1) == P_Source(index)*(dx/dt)+2*g_Pressure(index)- g_Pressure(index-1)-lowP+ lowP;
else
P_eqn(index) = 2*P_Corr(index) - P_Corr(index-1) - P_Corr(index+1)== P_Source(index)*(dx/dt)+ (2*g_Pressure(index))- g_Pressure(index-1)-g_Pressure(index+1);
end
end
[C,D] = equationsToMatrix(P_eqn,P_Corr);

% [C, D] = equationsToMatrix([P_eqn(1),P_eqn(2), P_eqn(3),P_eqn(4), P_eqn(5), P_eqn(6), P_eqn(7), P_eqn(8),P_eqn(9), P_eqn(10), P_eqn(11), P_eqn(12), P_eqn(13),P_eqn(14), P_eqn(15), P_eqn(16), P_eqn(17), P_eqn(18),P_eqn(19), P_eqn(20), P_eqn(21), P_eqn(22), P_eqn(23),P_eqn(24), P_eqn(25), P_eqn(26), P_eqn(27), P_eqn(28),P_eqn(29), P_eqn(30), P_eqn(31), P_eqn(32), P_eqn(33),P_eqn(34), P_eqn(35), P_eqn(36), P_eqn(37), P_eqn(38),P_eqn(39), P_eqn(40), P_eqn(41), P_eqn(42), P_eqn(43),P_eqn(44), P_eqn(45), P_eqn(46), P_eqn(47), P_eqn(48),P_eqn(49), P_eqn(50), P_eqn(51), P_eqn(52), P_eqn(53),P_eqn(54), P_eqn(55), P_eqn(56), P_eqn(57), P_eqn(58),P_eqn(59), P_eqn(60), P_eqn(61), P_eqn(62), P_eqn(63),P_eqn(64), P_eqn(65), P_eqn(66), P_eqn(67), P_eqn(68),P_eqn(69), P_eqn(70), P_eqn(71), P_eqn(72), P_eqn(73),P_eqn(74), P_eqn(75), P_eqn(76), P_eqn(77), P_eqn(78),P_eqn(79), P_eqn(80), P_eqn(81), P_eqn(82), P_eqn(83),P_eqn(84), P_eqn(85), P_eqn(86), P_eqn(87), P_eqn(88),P_eqn(89), P_eqn(90),P_eqn(91), P_eqn(92), P_eqn(93),P_eqn(94), P_eqn(95), P_eqn(96), P_eqn(97), P_eqn(98),P_eqn(99),P_eqn(99), P_eqn(100)], [ P_Corr(1), P_Corr(2), P_Corr(3), P_Corr(4), P_Corr(5), P_Corr(6), P_Corr(7), P_Corr(8), P_Corr(9), P_Corr(10), P_Corr(11), P_Corr(12), P_Corr(13), P_Corr(14), P_Corr(15), P_Corr(16), P_Corr(17), P_Corr(18), P_Corr(19), P_Corr(20), P_Corr(21), P_Corr(22), P_Corr(23), P_Corr(24), P_Corr(25), P_Corr(26), P_Corr(27), P_Corr(28), P_Corr(29), P_Corr(30), P_Corr(31), P_Corr(32), P_Corr(33), P_Corr(34), P_Corr(35), P_Corr(36), P_Corr(37), P_Corr(38), P_Corr(39), P_Corr(40), P_Corr(41), P_Corr(42), P_Corr(43), P_Corr(44), P_Corr(45), P_Corr(46), P_Corr(47), P_Corr(48), P_Corr(49), P_Corr(50), P_Corr(51), P_Corr(52), P_Corr(53), P_Corr(54), P_Corr(55), P_Corr(56), P_Corr(57), P_Corr(58), P_Corr(59), P_Corr(60), P_Corr(61), P_Corr(62), P_Corr(63), P_Corr(64), P_Corr(65), P_Corr(66), P_Corr(67), P_Corr(68), P_Corr(69), P_Corr(70), P_Corr(71), P_Corr(72), P_Corr(73), P_Corr(74), P_Corr(75), P_Corr(76), P_Corr(77), P_Corr(78), P_Corr(79), P_Corr(80), P_Corr(81), P_Corr(82), P_Corr(83), P_Corr(84), P_Corr(85), P_Corr(86), P_Corr(87), P_Corr(88), P_Corr(89), P_Corr(90), P_Corr(91), P_Corr(92), P_Corr(93), P_Corr(94), P_Corr(95), P_Corr(96), P_Corr(97), P_Corr(98), P_Corr(99), P_Corr(100)]);

P_1st_corr = linsolve(C,D);
% P_1st_corr(1) = 0;
% P_1st_corr(n_points-1) = 0;
% P_1st_Corr_new = g_Pressure + P_1st_corr';

% Corrected velocity
for index = 1: n_points
if index == 1
Corr_Vel(index) = 0;
elseif index == n_points
Corr_Vel(index) = 0;
else
% Corr_Vel(index) = Pred_Velo(index) + ((dt/(Den_on_Vnodes(index)*dx))*(P_1st_corr(index-1) - P_1st_corr(index)));
Corr_Vel(index) = ((dt/(Den_on_Vnodes(index)*dx))*(P_1st_corr(index-1) - P_1st_corr(index)));
end
end
% Corr_Vel_new = Corr_Vel + Pred_Velo;
% Mass flow
for index = 1:n_points-1
RHS_Mass(index) = ((DenS(index) -g_Density(index))*dx)/dt;
LHS_Mass(index) = (Den_on_Vnodes(index+1)*Corr_Vel(index+1)) - (Den_on_Vnodes(index)*Corr_Vel(index));
end

for index = 1:n_points
if index ==1
MachNumber(index) = 0;
elseif index ==n_points
MachNumber(index) = 0;
else
MachNumber(index) = Corr_Vel(index)./Sound_Speed(index);
end
end

guess_Vel = Corr_Vel;
g_Pressure = P_1st_corr;
g_Density= DenS;
time = time +1;
disp(time)
end

figure(2)
plot(x,Pred_Velo,'b','LineWidth',2)
hold on
plot(x,Corr_Vel,'r','LineWidth',2)
xlabel('x (m)')
ylabel('Velocity (m/s)')
legend('predicted velocity', 'Corrected velocity')
title ('Velocity')
saveas(gcf,'Velcity.png') % as matlab figure file to reopen in case
saveas(gcf,'Velocity.fig') % as matlab figure file to reopen in case


figure(4)
plot(x_Press, P_1st_corr ,'b','LineWidth',2)
hold on
plot(x_Press,Ini_P ,'r','LineWidth',2)
legend('Corrected P','Initial P')
xlabel('x (m)')
ylabel('Pressure (Pa)')
title ('Pressure')
saveas(gcf,'Pressure.png') % as matlab figure file to reopen in case
saveas(gcf,'Pressure.fig') % as matlab figure file to reopen in case