Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
Objective is to write a code to solve the Quasi 1D Subsonic-supersonic nozzle flow equations using the Macormack Method.
Initials Conditions:.
Where ,
Function for Non Conservative Form:
function [mfr_nc,P,Mach_nc,v,rho,T,Temp_throat,p_throat,v_throat,rho_throat, Mach_throat,mfr_throat]=non_conservative_form(x,dx,n,gamma,nt,c)
% Calculate intial profiles
rho=1-0.3146*x;
T=1-0.2341*x;
v=(0.1+1.09*x).*T.^0.5;
% Area
a=1+2.2*(x-1.5).^2;
throat=find(a==1);
% Determining the dt value
for m=1:n
dt_grid(m)=c.*(dx./(((T(m))^0.5+v(m))));
end
dt=min(dt_grid);
% Outer Loop
for k=1:nt
% Keeping an copy of initial values of rho,v,T for Comparision
rho_old=rho;
v_old=v;
T_old=T;
% Predictor Step
for j=2:n-1
% Spatial terms are discretized using FDS
dvdx=(v(j+1)-v(j))/dx;
drhodx=(rho(j+1)-rho(j))/dx;
dTdx=(T(j+1)-T(j))/dx;
dlog_a_dx=(log(a(j+1))-log(a(j)))/dx;
% Continuity Equation
drhodt_p(j)=-rho(j)*dvdx-rho(j)*v(j)*dlog_a_dx-v(j)*drhodx;
% Momentum Equation
dvdt_p(j)=-(v(j)*dvdx)-(1/gamma)*(dTdx+(T(j)*drhodx/rho(j)));
% Engergy Equation
dTdt_p(j)=-v(j)*dTdx-(gamma-1)*T(j)*(dvdx+v(j)*dlog_a_dx);
% 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
% Corretor Step
for j=2:n-1
% Spatial terms are discretized using BDS
dvdx=(v(j)-v(j-1))/dx;
drhodx=(rho(j)-rho(j-1))/dx;
dTdx=(T(j)-T(j-1))/dx;
dlog_a_dx=(log(a(j))-log(a(j-1)))/dx;
% Continuity Equation
drhodt_c(j)=-rho(j)*dvdx-rho(j)*v(j)*dlog_a_dx-(v(j)*drhodx);
% Momentum Equation
dvdt_c(j)=-v(j)*dvdx-(1/gamma)*(dTdx+(T(j)*drhodx/rho(j)));
% Engergy Equation
dTdt_c(j)=-v(j)*dTdx-(gamma-1)*T(j)*(dvdx+v(j)*dlog_a_dx);
end
% Compute the average time derivative
dvdt=0.5.*(dvdt_p+dvdt_c);
drhodt=0.5*(drhodt_p+drhodt_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 Contion
% Inlet BC
v(1)=2*v(2)-v(3);
% Outlet BC
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);
% Computing Pressure,Mass Flow Rate and Mach Number
P=rho.*T;
Mach_nc=(v./sqrt(T));
mfr_nc=rho.*a.*v;
% calculation of variables at the throat section of nozzle
rho_throat(k)=rho(throat);
Temp_throat(k)=T(throat);
v_throat(k)=v(throat);
p_throat(k)=P(throat);
Mach_throat(k)=Mach_nc(throat);
mfr_throat(k)=mfr_nc(throat);
% Plotting mass flow rate at different timesteps
figure(1)
if k==50
plot(x,mfr_nc,'b','linewidth',1);
hold on
elseif k==100
plot(x,mfr_nc,'r','linewidth',1);
hold on
elseif k==150
plot(x,mfr_nc,'c','linewidth',1);
hold on
elseif k==200
plot(x,mfr_nc,'g','linewidth',1);
hold on
elseif k==300
plot(x,mfr_nc,'y','linewidth',1);
elseif k==400
plot(x,mfr_nc,'k','linewidth',1);
elseif k==700
plot(x,mfr_nc,'*','linewidth',1);
elseif k==1400
plot(x,mfr_nc,'b','linewidth',1);
grid on
xlabel('Non Dimesional Non Conservative Distance ','fontweight','bold','fontsize',8);
ylabel('Non Dimesional Non Conservative Mass Flow Rate','fontweight','bold','fontsize',8);
title('Non Dimesional Non Conservative Mass Flow Rate at different timesteps','fontsize',8,'color','r');
legend({'50th timestep','100th timestep','150th timestep','200th timestep','300th timestep','400th timestep','700th timestep','1400th timestep'});
end
end
end
Function for Conservative Form:
function [mfr_c,P_c,Mach_c,v_c,rho_c,T_c,Temp_throat_c,p_throat_c,v_throat_c,rho_throat_c, Mach_throat_cs,mfr_throat_c]=conservative_form(x,dx,c,n,gamma,nt)
rho_c=ones(1,n);
T_c=ones(1,n);
% Area
a=1+2.2*(x-1.5).^2;
% Calculate intial profiles
for i=1:n
if (x(i)>=0 && x(i)<=0.5)
rho_c(i)=1;
T_c(i)=1;
elseif (x(i)>=0.5 && x(i)<=1.5)
rho_c(i)=1-0.366*(x(i)-0.5);
T_c(i)=1-0.167*(x(i)-0.5);
elseif (x(i)>=1.5 && x(i)<=3)
rho_c(i)=0.634-0.3879*(x(i)-1.5);
T_c(i)=0.833-0.3507*(x(i)-1.5);
end
end
P_c=rho_c.*T_c;
v_c=0.59./(rho_c.*a);
c=0.5;
for i=1:n
dt_grid(i)=c.*(dx./((((T_c(i).^0.5))+v_c(i))));
end
dt=min(dt_grid);
% dt=0.0267;
% Initial Conditions
U1=rho_c.*a;
U2=U1.*v_c;
U3=U1.*((T_c./(gamma-1))+(gamma/2).*v_c.^2);
throat=find(a==1);
% Outer Loop
for k=1:nt
U1_old=U1;
U2_old=U2;
U3_old=U3;
% Determining Flux terms for predictor step
F1=U2;
F2=((U2.^2)./U1)+((gamma-1)/gamma)*(U3-(gamma/2)*((U2.^2)./U1));
F3=((gamma*U2.*U3)./U1)-(gamma*(gamma-1)*((U2.^3)./(2*U1.^2)));
% Predictor Step (FDS)
for l=2:n-1
% dadx(l)=((a(l+1))-(a(l)))/dx;
% j2(l)=((gamma-1)/gamma).*(U3(l)-(gamma/2)*(U2(l).^2)./U1(l))*dlogadx(l);
j2(l)=(1/gamma)*rho_c(l)*T_c(l)*((a(l+1)-a(l))/dx);
% Continuity Equation
dU1dt_p(l)=-((F1(l+1)-F1(l))/dx);
% Momentum Equation
dU2dt_p(l)=-((F2(l+1)-F2(l))/dx)+j2(l);
% Engergy Equation
dU3dt_p(l)=-((F3(l+1)-F3(l))/dx);
% Solution Update
U1(l)=U1(l)+dU1dt_p(l)*dt;
U2(l)=U2(l)+dU2dt_p(l)*dt;
U3(l)=U3(l)+dU3dt_p(l)*dt;
end
% flux update
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 BDS
for l=2:n-1
% dadx(l)=((a(l)-a(l-1)))/dx;
% dlogadx(l)=(log(a(l))-(log(a(l-1))))/dx;
% j2(l)=((gamma-1)/gamma).*(U3(l)-(gamma/2)*(U2(l).^2)./U1(l))*dlogadx(l);
j2(l)=(1/gamma)*rho_c(l)*T_c(l)*((a(l)-a(l-1))/dx);
% Continuity Equation
dU1dt_c(l)=-((F1(l)-F1(l-1))/dx);
% Momentum Equation
dU2dt_c(l)=-((F2(l)-F2(l-1))/dx)+j2(l);
% Engergy Equation
dU3dt_c(l)=-((F3(l)-F3(l-1))/dx);
end
% Average Time Derivative
dU1dt=0.5*(dU1dt_p+dU1dt_c);
dU2dt=0.5*(dU2dt_p+dU2dt_c);
dU3dt=0.5*(dU3dt_p+dU3dt_c);
% Final Solution Update
for j=2:n-1
U1(j)=U1_old(j)+dU1dt(j)*dt;
U2(j)=U2_old(j)+dU2dt(j)*dt;
U3(j)=U3_old(j)+dU3dt(j)*dt;
end
% Boundary Conditons
% Inlet
U1(1)=rho_c(1)*a(1);
U2(1)=2*U2(2)-U2(3);
v_c(1)=U2(1)./U1(1);
U3(1)=U1(1)*((T_c(1)/(gamma-1))+((gamma/2)*(v_c(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);
% Non Dimesional Properties
rho_c=U1./a;
v_c=U2./U1;
T_c=(gamma-1)*((U3./U1- (gamma/2)*(v_c).^2));
% Computing Pressure,Mass Flow Rate and Mach Number
P_c=rho_c.*T_c;
Mach_c=(v_c./sqrt(T_c));
mfr_c=rho_c.*v_c.*a;
% calculation of variables at the throat section of nozzle
rho_throat_c(k)=rho_c(throat);
Temp_throat_c(k)=T_c(throat);
v_throat_c(k)=v_c(throat);
p_throat_c(k)=P_c(throat);
Mach_throat_cs(k)=Mach_c(throat);
mfr_throat_c(k)=mfr_c(throat);
% Plotting mass flow rate at different timesteps
figure(2)
if k==50
plot(x,mfr_c,'b','linewidth',1);
hold on
elseif k==100
plot(x,mfr_c,'r','linewidth',1);
hold on
elseif k==150
plot(x,mfr_c,'c','linewidth',1);
hold on
elseif k==200
plot(x,mfr_c,'g','linewidth',1);
hold on
elseif k==300
plot(x,mfr_c,'y','linewidth',1);
elseif k==400
plot(x,mfr_c,'k','linewidth',1);
elseif k==700
plot(x,mfr_c,'*','linewidth',1);
elseif k==1400
plot(x,mfr_c,'b','linewidth',1);
grid on
xlabel('Non Dimesional Conservative Distance ','fontweight','bold','fontsize',8);
ylabel('Non Dimesional Conservative Mass Flow Rate','fontweight','bold','fontsize',8);
title('Non Dimesional Conservative Mass Flow Rate at different timesteps','fontsize',8,'color','r');
legend('50th timestep','100th timestep','150th timestep','200th timestep','300th timestep','400th timestep','700th timestep','1400th timestep');
end
end
Main Run Code:
close all
clear all
clc
% Inputs
n=31;
x=linspace(0,3,n);
dx=x(2)-x(1);
gamma=1.4;
% Throat Point
% No. of timestep
nt=1400;
c=0.5; % Courant number
tic;
% Function For Non Conservative Form
[mfr_nc,P,Mach_nc,v,rho,T,Temp_throat,p_throat,v_throat,rho_throat, Mach_throat,mfr_throat]=non_conservative_form(x,dx,n,gamma,nt,c);
time_elapsed_non_cons=toc;
fprintf('The Computational Time for non conservative form is: %0.4f',time_elapsed_non_cons);
hold on
tic;
% Function For Non Conservative Form
[mfr_c,P_c,Mach_c,v_c,rho_c,T_c,Temp_throat_c,p_throat_c,v_throat_c,rho_throat_c, Mach_throat_cs,mfr_throat_c]=conservative_form(x,dx,c,n,gamma,nt);
time_elapsed_cons=toc;
fprintf('The Computational Time for conservative form is: %0.4f s',time_elapsed_cons);
% Plotting Information along throat
figure(3)
a1=subplot(4,1,1);
plot(linspace(1,nt,nt),rho_throat,'color','r')
hold on
plot(linspace(1,nt,nt),rho_throat_c,'color','b')
xlabel('Timestep');
ylabel('Density ratio at throat');
legend(a1,{'Non Conservative Form','Conservative Form'});
title('Time-wise Variation of flow field variables at throat');
grid minor
a2=subplot(4,1,2);
plot(linspace(1,nt,nt),Temp_throat,'color','r')
hold on
plot(linspace(1,nt,nt),Temp_throat_c,'color','b')
xlabel('Timestep');
ylabel('Temperature ratio at throat');
legend(a2,{'Non Conservative Form','Conservative Form'});
grid minor
a3=subplot(4,1,3);
plot(linspace(1,nt,nt),p_throat,'color','r')
hold on
plot(linspace(1,nt,nt),p_throat_c,'color','b')
xlabel('Timestep');
ylabel('Pressure ratio at throat');
legend(a3,{'Non Conservative Form','Conservative Form'});
grid minor
a4=subplot(4,1,4);
plot(linspace(1,nt,nt),Mach_throat,'color','r')
hold on
plot(linspace(1,nt,nt),Mach_throat_cs,'color','b')
xlabel('Timestep');
ylabel('Mach Number at throat');
legend(a4,{'Non Conservative Form','Conservative Form'});
grid minor
% Comparision of variation of flow field variables along x axis
figure(4)
b1=subplot(4,1,1);
plot(x,rho,'color','g')
hold on
plot(x,rho_c,'*')
xlabel('Non Dimensional nozzle length');
ylabel('Density ratio');
legend(b1,{'Non Conservative Form','Conservative Form'});
title('Variation of flow field variables along x axis');
grid minor
b2=subplot(4,1,2);
plot(x,T,'color','g')
hold on
plot(x,T_c,'*')
xlabel('Non Dimensional nozzle length');
ylabel('Temperature ratio ');
legend(b2,{'Non Conservative Form','Conservative Form'});
grid minor
b3=subplot(4,1,3);
plot(x,P,'color','g')
hold on
plot(x,P_c,'*')
xlabel('Non Dimensional nozzle length');
ylabel('Pressure Ratio ');
legend(b3,{'Non Conservative Form','Conservative Form'});
grid minor
hold on
b4=subplot(4,1,4);
plot(x,Mach_nc,'color','g')
hold on
plot(x,Mach_c,'*')
xlabel('Non Dimensional nozzle length');
ylabel('Mach Number ');
legend(b4,{'Non Conservative Form','Conservative Form'});
grid minor
hold on
% Comparision of mass flow rate between conservative form and non-conservative form
figure(5);
hold on
plot(x,mfr_c,'g','linewidth',1);
hold on;
plot(x,mfr_nc,'y','linewidth',1);
hold on;
legend('Conservative Form','Non Conservative Form');
xlabel('Non Dimensional Length of the Domain');
ylabel('Non Dimensionbal Mass flow rate ');
title('Comparision of Mass flow rates')
grid minor
Steady-state distribution of primitive variables inside the nozzle:
Time-wise variation of the primitive variables:
Variation of Mass flow rate distribution inside the nozzle at different time steps during the time-marching process:
For Non Conservative form:
For Conservative form:
Comparison of Normalized mass flow rate distributions of both forms: