Week 7 - Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method

clear all
close all
clc

% inputs
n= 31
x= linspace(0,3,n)
dx= x(2)-x(1);
gamma = 1.4;

% calculate initail profiles
rho = 1-0.3146*x;
t=1-0.2314*x; % t= temperature
v=(0.1+1.09*x).*t.^0.5;

% area
a= 1+2.2*(x-1.5).^2;

% time steps
nt = 1400;
dt=0.01;

% outer time loop
for k= 1:nt
rho_old=rho;
v_old=v;
t_old = t;

% predictor method
for j= 2:n-1
%drho_dt_p=-rho(j)*(v(j+1)-v(j))/dx-rho(j)*v(j)*(log(a(j+1))-log(a(j)))/dx-v(j)*(rho(j+1)-rho(j))/dx;
dvdx= (v(j+1)-v(j))/dx;
drhodx = (rho(j+1)-rho(j))/dx;
dtdx = (t(j+1)-t(j))/dx;
dlogadx = (log(a(j+1))-log(a(j)))/dx;

% continutity equation

drhodt_p(j) = -rho(j)*dvdx - rho(j)*v(j)*dlogadx-v(j)*drhodx;
% momentum equation
dvdt_p(j)= - v(j)*dvdx -(1/gamma)*(dtdx +(t(j)/rho(j))*drhodx);


% energy equation
dtdt_p(j)= -v(j)*dtdx -(gamma-1)*(dvdx + v(j)*dlogadx);


% solution update
rho(j)= rho(j)+drhodt_p(j)*dt;
v(j)= v(j) + dvdt_p(j)*dt;
t(j)=t(j)+dtdt_p(j)*dt;
end
%corrector method
for j= 2:n-1
%drho_dt_p=-rho(j)*(v(j+1)-v(j))/dx-rho(j)*v(j)*(log(a(j+1))-log(a(j)))/dx-v(j)*(rho(j+1)-rho(j))/dx;
dvdx= (v(j)-v(j-1))/dx;
drhodx = (rho(j)-rho(j-1))/dx;
dtdx = (t(j)-t(j-1))/dx;
dlogadx = (log(a(j))-log(a(j-1)))/dx;

% continutity equation

drhodt_c(j) = -rho(j)*dvdx - rho(j)*v(j)*dlogadx-v(j)*drhodx;
% momentum equation
dvdt_c(j)= - v(j)*dvdx -(1/gamma)*(dtdx +(t(j)/rho(j))*drhodx);


% energy equation
dtdt_c(j)= -v(j)*dtdx -(gamma-1)*(dvdx + v(j)*dlogadx);


end
% compute the average time derivative
drhodt= 0.5* (drhodt_p + drhodt_c);
dvdt = 0.5*(dvdt_p + dvdt_c);
dtdt = 0.5*(dtdt_p + dtdt_c);

% final solution update
for i = 2:n-1
rho(i) = rho_old(i)+ drhodt(i)*dt;
v(i)= v_old(i)+ dvdt(i)*dt;
t(i)= t_old(i)+ dtdt(i)*dt;
end

%Apply boundary conditions
% inlet
v(1)=2*v(2)- v(3);
%outlet
v(n) = 2*v(n-1)-v(n-2);
rho(n)=2*rho(n-1)-rho(n-2);
t(n)=2*t(n-1)-t(n-2);

p = rho.*t;
% plot(x,v)
% title(sprintf('time step=%d',k))
% pause(0.03)
end

mach=(v)./(t.^0.5);

figure(1)
subplot(3,1,1)
plot(x,rho)
ylabel('density_ratio')
subplot(3,1,2)
plot(x,t)
ylabel('temperature_ratio')
subplot(3,1,3)
plot(x,p)
ylabel('pressure')