Solution of quasi static1D supersonic nozzle flow equations conservative and non conservative using Macromack method for

Objective: To solve the conservative and non conservative governing equations of 1D supersonic nozzle and understand its flow properties at throat,and to perform a grid independence test

Code:

Code for Non conservative form

close all
clear all
clc
n=31;
gamma=1.4;
x =linspace(0,3,n);
dx=x(2)-x(1);
a=1+2.2*(x-1.5).^(2);
th=find(a==1);
rho=1-0.3146*x;
T = 1-0.231*x;
v=(0.1+1.09*x).*T.^(1/2);
p=rho.*T;
%plot(x,a)
c=0.5
nt=1000;
for k=1:nt;
dt = min(c*(dx./(sqrt(T)+v)));
rho_old =rho;
v_old=v;
T_old=T;
% predctor method
for i=2:n-1;
% continuity equation
dv_dx = (v(i+1)-v(i))/dx;
dlog_a_dx =(log(a(i+1))-log(a(i)))/dx;
drho_dx = (rho(i+1)-rho(i))/dx;
dT_dx= (T(i+1)-T(i)) /dx;
drho_dt_p(i) = -rho(i)*(dv_dx) - rho(i)*v(i)*(dlog_a_dx) - v(i)*(drho_dx);
dv_dt_p(i) = -v(i)*(dv_dx)-(1/gamma)*((dT_dx)+(T(i)/rho(i))*drho_dx);
dT_dt_p(i) = -v(i)*(dT_dx) - (gamma-1)*T(i)*(dv_dx +v(i)*dlog_a_dx);
rho(i)=rho(i)+drho_dt_p(i)*(dt);
v(i)=v(i)+dv_dt_p(i)*(dt);
T(i)=T(i)+dT_dt_p(i)*(dt);
end
%CORRECTOR METHOD
for j =2:n-1
dv_dx = (v(j)-v(j-1))/dx;
dlog_a_dx =(log(a(j))-log(a(j-1)))/dx;
drho_dx = (rho(j)-rho(j-1))/dx;
dT_dx= (T(j)-T(j-1)) /dx;
drho_dt_c(j) = -rho(j).*(dv_dx) - rho(j)*v(j)*(dlog_a_dx) - v(j)*(drho_dx);
dv_dt_c(j) = -v(j)*dv_dx - (1/gamma)*(dT_dx + (T(j)/rho(j))*drho_dx);
dT_dt_c(j) = -v(j)*dT_dx - (gamma-1)*T(j)*(dv_dx +v(j)*dlog_a_dx);
end
drho_dt=0.5*(drho_dt_c + drho_dt_p);
dv_dt=0.5*(dv_dt_c + dv_dt_p);
dT_dt=0.5*(dT_dt_c + dT_dt_p);
for i=2:n-1
rho(i)=rho_old(i)+drho_dt(i)*dt;
v(i)=v_old(i)+dv_dt(i)*dt;
T(i)=T_old(i)+dT_dt(i)*dt;
end
%boundary condition
%inlet
v(1) =2*v(2)-v(3);
%outlet
rho(n) =2*rho(n-1)-rho(n-2);
v(n) =2*v(n-1)-v(n-2);
T(n) =2*T(n-1)-T(n-2);
M =v./ sqrt(T);
p=rho.*T;
mass =rho.*v.*a;
M_throat(k) = M(th);
p_throat(k) = p(th);
T_throat(k) = T(th);
rho_throat(k) = rho(th);
plot_nonconservative_diferent_time(x,k,mass)
end
% plot_nonconservative_diferent_time(x,k,mass)
figure(1)
hold on
plot(T_throat,\'b\')
plot(p_throat,\'g\')
plot(rho_throat,\'magenta\')
plot(M_throat,\'black\')
legend(\'T throat\',\'p throat\',\'rho throat\',\'M throat\')
xlabel(\'no. of iteration\');
title(sprintf(\'Variaation of properties at throat section \\n NON CORSERVATIVE FORM\'))
figure(2)
subplot(4,1,1)
plot(x,M)
ylabel(\'M--->\')
title(sprintf(\'ONE DIMENSIONAL SUPERSONIC NOZZLE FLOW \\n (NON CONSERVATION FORM) \\n
NO. OF TIME STEP=%d\',k))
subplot(4,1,2)
plot(x,T)
% xlabel(\'LENGTH--->\')
ylabel(\'Temprature ratio--->\')
subplot(4,1,3)
plot(x,rho)
% xlabel(\'NON DIMENSIONAL LENGTH-->\')
ylabel(\'Density ratio-->\')
subplot(4,1,4)
plot(x,p)
xlabel(\'NON DIMENSIONAL LENGTH-->\')
ylabel(\'pressure ratio-->\')
pause(0.03)

Code for plotting non conservative result at different time steps

function[] = plot_nonconservative_diferent_time(x,k,mass)
figure(3)
grid on
if k==1
plot(x,mass,\'magenta\')
hold on
elseif k==100
plot(x,mass,\'b\')
elseif(k==300)
plot(x,mass,\'cyan\')
elseif(k==500)
plot(x,mass,\'g\')
elseif(k==700)
plot(x,mass,\'y\')
elseif(k==850)
plot(x,mass,\'r\')
elseif(k==1000)
plot(x,mass,\'black\')
legend(\'k=1\',\'k=100\',\'k=300\',\'k=500\',\'k=700\',\'k=850\',\'k=1000\')
xlabel(\'non dimensional length\')
ylabel(\'non dimensional mass flow rate\')
title(\'Variation in mass flow rate at diffetent time step ( NON CONSERVATIVE FORM)\')
end
end

Code for conservative form

close all
clear all
clc
n=31;
x=linspace(0,3,n);
dx=x(2)-x(1);
a= 1+2.2*(x-1.5).^(2);
th= find(a==1);
for i=1:n
if (0<= x(i))&&(x(i)<= 0.5)
rho(i)=1;
T(i)=1;
elseif (0.5 < x(i))&&(x(i)<= 1.5)
rho(i)=1-(0.366*(x(i)-0.5));
T(i) = 1-0.167*(x(i)-0.5);
else
rho(i) =0.634-(0.3879*(x(i)-1.5));
T(i) = 0.833-(0.3507*(x(i)-1.5));
end
end
u1=rho.*a;
% u2=0.590
gamma=1.4;
v=0.590./(u1);
u2=u1.*v;
u3=u1.*((T/(gamma-1)+ (0.5*gamma*v.^2)));
c=0.5;
for t=1:1500
dt = min(c*(dx./(sqrt(T)+v)));
u1_old=u1;
u2_old=u2;
u3_old=u3;
%Flux variables for predictor step
f1 = u2;
f2 = ((u2.^2)./u1) + (((gamma-1)/gamma)*(u3 - ((gamma*0.5*(u2.^2))./u1)));
f3 = ((gamma*u2.*u3)./u1) - ((gamma*0.5*(gamma-1)*(u2.^3))./(u1.^2));
%Predictor step
for i=2:n-1
a_dx=(a(i+1)-a(i))/dx;
j2(i)=(rho(i).*T(i).*a_dx)/gamma;
du1_dt_p(i)= -(f1(i+1)-f1(i))/dx;
du2_dt_p(i)= -(f2(i+1)-f2(i))/dx + j2(i);
du3_dt_p(i) = -(f3(i+1)-f3(i))/dx;
u1(i) = u1_old(i) + du1_dt_p(i)*dt;
u2(i) = u2_old(i) + du2_dt_p(i)*dt;
u3(i) = u3_old(i) + du3_dt_p(i)*dt;
%find premitive variable for corrector term
rho(i)=u1(i)./a(i);
v(i)=u2(i)./u1(i);
T(i) = (gamma-1)*((u3(i)./u1(i))-(0.5*gamma)*(v(i)^2));
end
%Updated Flux variables for corrector step
f1 = u2;
f2 = ((u2.^2)./u1) + (((gamma-1)/gamma)*(u3 - ((gamma*0.5*(u2.^2))./u1)));
f3 = ((gamma*u2.*u3)./u1) - ((gamma*0.5*(gamma-1)*(u2.^3))./(u1.^2));
%Corrector step
for i=2:n-1
a_dx=(a(i)-a(i-1))/dx;
j2(i)=(rho(i).*T(i).*a_dx)/gamma;
du1_dt_c(i)= -(f1(i)-f1(i-1))/dx;
du2_dt_c(i)= -(f2(i)-f2(i-1))/dx + j2(i);
du3_dt_c(i) = -(f3(i)-f3(i-1))/dx;
end
% Average values
du1_dt =0.5*((du1_dt_p) + (du1_dt_c));
du2_dt =0.5*((du2_dt_p) + (du2_dt_c)) ;
du3_dt =0.5*((du3_dt_p) + (du3_dt_c));
for i =2:n-1
u1(i)= u1_old(i) + du1_dt(i)*dt;
u2(i)= u2_old(i) + du2_dt(i)*dt;
u3(i)= u3_old(i) + du3_dt(i)*dt;
end
%Final update of primitive variables
rho = u1./a;
v = u2./u1;
T = (gamma-1)*((u3./u1)-((0.5*gamma)*(v.^2)));
%Boundary conditions
%Inlet
u1(1) = rho(1)*a(1);
u2(1) = 2*u2(2)-u2(3);
u3(1) = u1(1)*((T(1)/(gamma-1))+(0.5*gamma*(v(1).^2)));
%Outlet
u1(n) = 2*u1(n-1)-u1(n-2);
u2(n) = 2*u2(n-1)-u2(n-2);
u3(n) = 2*u3(n-1)-u3(n-2);
p = rho.*T;
mass= rho.*v.*a;
M = v./(sqrt(T));
%CALCUTING PROPERTIES AT THROAT SECTION
T_throat(t) = T(th);
M_throat(t)=M(th);
rho_throat(t)=rho(th);
v_throat=v(th);
mass_throat = rho_throat.*v_throat.*a(th);
p_throat=rho_throat.*T_throat;
plot_conservative_diff_time(mass,x,t)
figure(1)
end
plot_conservative_throat(T_throat,M_throat,rho_throat,p_throat)
plot_conservation(x,M,T,rho,p,t)

Code for plotting conservative solution at different time

function[] = plot_conservative_diff_time(mass,x,t)
figure(1)
if t==1
plot(x,mass,\'r\')
hold on
elseif t==100
plot(x,mass,\'b\')
elseif(t==300)
plot(x,mass)
elseif(t==500)
plot(x,mass,\'g\')
elseif(t==700)
plot(x,mass,\'y\')
elseif(t==1000)
plot(x,mass,\'black\')
elseif(t==1400)
plot(x,mass,\'magenta\')
legend(\'t=1\',\'t=100\',\'t=300\',\'t=500\',\'t=700\',\'t=1000\',\'t=1400\')
xlabel(\'non dimensional length\')
ylabel(\'non dimensional mass flow rate\')
title(\'Variation in mass flow rate at diffetent time step ( CONSERVATIVE FORM)\')
end
end

Code for plotting conservative results at throat

function[] = plot_conservative_throat(T_throat,M_throat,rho_throat,p_throat)
figure(2)
plot(T_throat,\'g\')
hold on
plot(M_throat,\'r\')
plot(rho_throat,\'b\')
plot(p_throat,\'m\')
legend(\'T throat\',\'M throat\',\'rho throat\',\'p throat\')
xlabel(\'no. of iteration\')
title(sprintf(\'Variation of properties according to iteration at throat position \\n
CONSERVATIVE FORM\'))
end

Code for conservation function:

function[] = plot_conservation(x,M,T,rho,p,t)
grid on
figure(3)
subplot(4,1,1)
plot(x,M)
ylabel(\'M--->\')
title(sprintf(\'ONE DIMENSIONAL SUPERSONIC NOZZLE FLOW \\n ( CONSERVATION FORM) \\n NO. OF T
IME STEP=%d\',t))
subplot(4,1,2)
plot(x,T)
% xlabel(\'LENGTH--->\')
ylabel(\'Temprature ratio--->\')
subplot(4,1,3)
plot(x,rho)
% xlabel(\'NON DIMENSIONAL LENGTH-->\')
ylabel(\'Density ratio-->\')
subplot(4,1,4)
plot(x,p)
%xlabel(\'NON DIMENSIONAL LENGTH-->\')
ylabel(\'pressure ratio-->\')
pause(0.03)
end

Results:

Non conservative results

Grid independence test results for N= 31 and N=62 respectively

For N=31

For N=62

Conservative form results

Grid indepence test results for conservative form for N=31 and N=62 respectively

For N=31

For N=62

1) For N= 62 Non conservative solution converges between 1300 and 1400 iterations while conservative  solution takes more tham 1400 iterations to converge. For N=31, non conservative solution converges between 700  and 850 iteration while in conservative solution takes between 1000 and 1100

2) Grid independece is achieved at 124 no of grid points after which solution does not change

3) For same no of grid points, conservative solution gives result closer to analytical solution compared to non conservative solution