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

Aim:-

Simulation of a 1D Super-sonic nozzle flow using Macormack Method for
conservative form as well as for non-conservative form.
Discretization:-The process of converting the PDE equation into algebric equation is
basically discretization.Now, here in this step, we will use the Mccormack method to solve it.
Boundary conditions:-

Inlet:-
u1(1)=ρ(1)*a(1);
u2(1)=2*u2(2)-u2(3);
u3(1)=u1(1)*((t(1)/(γ-1))+((γ/2)*v(1)^2));
v(1)=u2(1)/u1(1);
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);
function [rho_con,v_con,t_con,mfl_con,m_con,p_con]=Conservation(n,nt)
%creating mesh & assiging input
x=linspace(0,3,n);
gamma=1.4;
dx=x(2)-x(1);
c =0.5;
%calculate the initial profiles
for i=1:n
if (x(i)>=0 && x(i)<=0.5)
rho_con(i)=1
t_con(i)=1;
elseif (x(i)>0.5 && x(i)<1.5)
rho_con(i)=1-0.366.*(x(i)-0.5);
t_con(i)=1.0-0.167*(x(i)-0.5);
elseif (x(i)>=1.5 && x(i)<=3.5)
rho_con(i)=0.634-0.3879*(x(i)-1.5)
t_con(i)=0.833-0.3507*(x(i)-1.5)
end
end
a=1+2.2*(x-1.5).^2; %Area
v_con=0.59./(rho_con.*a); %velocity
%Defining solution vector
u1=rho_con.*a;
u2=rho_con.*a.*v_con;
u3= rho_con.*a.*((t_con/(gamma-1))+((gamma/2)*v_con.^2));
f1 = u2;
f2=(((u2).^2)./(u1)) +
((gamma-1)/gamma)*(u3-((gamma/2)*((u2.^2)./(u1))));
f3=((gamma*u2.*u3)./(u1))-((gamma*(gamma-1))/2).*((u2.^3)./(u1.^2));
%outer loop
for k=1:nt
u1_old=u1;
u2_old=u2;
u3_old=u3;
dt = min(c*(dx./(sqrt(t_con)+v_con)));
% predicator method
for j=2:n-1
f1(j) = u2(j);
f2(j)=(((u2(j))^2)/(u1(j))) +
((gamma-1)/gamma)*(u3(j)-((gamma/2)*((u2(j)^2)/(u1(j)))));
f3(j)=((gamma*u2(j)*u3(j))/(u1(j)))-((gamma*(gamma-1))/2)*((u2(j)^3)/(u1(j
)^2));
df1dx(j)=(f1(j+1)-f1(j))/dx;
df2dx(j)=(f2(j+1)-f2(j))/dx;
df3dx(j)=(f3(j+1)-f3(j))/dx;
dlogadx(j)=((a(j+1)-a(j)))/dx;
J2=(1/gamma)*rho_con(j).*t_con(j).*dlogadx(j);
% continuity equation
du1dt_p(j)=-df1dx(j);
%Momentum equation
du2dt_p(j)=-df2dx(j)+J2;
% energy equation
du3dt_p(j)=-df3dx(j);
% solution update
u1(j)=u1(j)+du1dt_p(j)*dt;
u2(j)=u2(j)+du2dt_p(j)*dt;
u3(j)=u3(j)+du3dt_p(j)*dt;
rho_con(j)=u1(j)/a(j);
v_con(j)=u2(j)/u1(j);
t_con(j)=(gamma-1)*((u3(j)/u1(j))-(0.5*gamma*v_con(j)^2));
end
f1 = u2;
f2=(((u2).^2)./(u1)) +
((gamma-1)/gamma).*(u3-((gamma/2)*((u2.^2)./(u1))));
f3=((gamma*u2.*u3)./(u1))-((gamma*(gamma-1))/2).*((u2.^3)./(u1.^2));
% corrector method
for j=2:n-1
df1dx(j)=(f1(j)-f1(j-1))/dx;
df2dx(j)=(f2(j)-f2(j-1))/dx;
df3dx(j)=(f3(j)-f3(j-1))/dx;
dlogadx(j)=((a(j)-a(j-1)))/dx;
J2=(1/gamma)*rho_con(j).*t_con(j).*dlogadx(j);
% continuity equation
du1dt_c(j)=-df1dx(j);
%Momentum equation
du2dt_c(j)=-df2dx(j) + J2;
% energy equation
du3dt_c(j)=-df3dx(j);
end
% compute the average period
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 i=2:n-1
u1(i)=u1_old(i)+du1dt(i)*dt;
u2(i)=u2_old(i)+du2dt(i)*dt;
u3(i)=u3_old(i)+du3dt(i)*dt;
end
% calculating the flow parameters
rho_con=u1./a;
v_con=u2./u1;
t_con=(gamma-1)*((u3./u1)-(0.5.*gamma*v_con.^2));
% apply boundary conditions
% inlet boundary
u1(1) = rho_con(1)*a(1);
u2(1) = 2*u2(2) - u2(3);
u3(1) = u1(1)*((t_con(1)/(gamma-1))+(0.5*gamma*(v_con(1)^2)));
% outlet boundary
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);
% calculation of the variables
m_con=v_con./sqrt(t_con); % Mach number
mfl_con=rho_con.*v_con.*a; % mass flow
p_con=rho_con.*t_con; % pressure
%Plotting mass flow rate at different time steps
figure(2)
if k==1
plot(x,mfl_con,'*','color','b','markersize',10)
hold on
elseif k==100
plot(x,mfl_con,'y','linewidth',1.5)
hold on
elseif k==200
plot(x,mfl_con,'r','linewidth',1.5)
hold on
elseif k==700
plot(x,mfl_con,'c','linewidth',1.5)
hold on
elseif k==1000
plot(x,mfl_con,'--','color','b','linewidth',1.5)
hold on
xlabel('Non_dimensional distance x')
ylabel('Mass flow rate')
legend('1 time step','100 time steps','200 time steps','700 time
steps','1000 time steps')
legend('Location','north')
title(['Conversation analysis';'Mass flow rate at different time
steps'])
grid on
end
end
end

#Non-Conservation form:- Non-Conservation form refers to an arrangement of an
equation or system of equations, usually representing a hyperbolic, that emphasizes
that a property represented is non-conserved.
Initial profile:-In this step we assign initial profile for different parameters
● ρ=1-0.3146*x;
● v=(0.1+1.09*x).*t.^0.5;
● a=1+2.2*(x-1.5).^2;
● t=1-0.2314*x;
Discretization:-The process of converting the PDE equation into algebric equation is
basically discretization. Now, here in this step, we will use the Mccormack method to solve it.
Boundary conditions:-
Inlet:-
v(1)=2*v(2)-v(3);
t(1)=2*t(2)-t(3);
ρ(1)=2*ρ(2)-ρ(3);
Outlet:-
v(n)=2*v(n-1)-v(n-2);
t(n)=2*t(n-1)-t(n-2);
ρ(n)=2*ρ(n-1)-ρ(n-2);
function [rho_nc,v_nc,t_nc,mfl_nc,M_nc,p_nc]=Non_Conservation(n,nt)
%Solving 1D qaussi supersonic equation using maccormak method for
conservative form
%Creating mesh and assigning inputs
x=linspace(0,3,n);
dx=x(2)-x(1);
gamma=1.4;
%Initial profile
rho_nc=1-0.3146*x;
t_nc=1-0.2314*x;
v_nc=(0.1+1.09*x).*t_nc.^0.5;
%Area
a=1+2.2*(x-1.5).^2;
%time steps
dt=1e-3;%time steps size
%time loop
for k=1:nt
%Old values
rho_nc_old=rho_nc;
t_nc_old=t_nc;
v_nc_old=v_nc;
%Predicitor step
for j=2:n-1 %space loop
drhodx=(rho_nc(j+1)-rho_nc(j))/dx;
dvdx=(v_nc(j+1)-v_nc(j))/dx;
dtdx=(t_nc(j+1)-t_nc(j))/dx;
dlogadx=(log(a(j+1))-log(a(j)))/dx;
%Continuity equation
drhodt_p(j)=-rho_nc(j)*dvdx-rho_nc(j)*v_nc(j)*dlogadx-v_nc(j)*drhodx;
%Energy equation
dtdt_p(j)=-v_nc(j)*dtdx-(gamma-1)*t_nc(j)*(dvdx+v_nc(j)*dlogadx);
%Momentum equation
dvdt_p(j)=-v_nc(j)*dvdx-(1/gamma)*(dtdx+(t_nc(j)/rho_nc(j))*drhodx);
%Solution update
v_nc(j)=v_nc(j)+dvdt_p(j)*dt;
rho_nc(j)=rho_nc(j)+drhodt_p(j)*dt;
t_nc(j)=t_nc(j)+dtdt_p(j)*dt;
end
%Corrector step
for j=2:n-1 %space loop
drhodx=(rho_nc(j)-rho_nc(j-1))/dx;
dvdx=(v_nc(j)-v_nc(j-1))/dx;
dtdx=(t_nc(j)-t_nc(j-1))/dx;
dlogadx=(log(a(j))-log(a(j-1)))/dx;
%Continuity equation
drhodt_c(j)=-rho_nc(j)*dvdx-rho_nc(j)*v_nc(j)*dlogadx-v_nc(j)*drhodx;
%Energy equation
dtdt_c(j)=-v_nc(j)*dtdx-(gamma-1)*t_nc(j)*(dvdx+v_nc(j)*dlogadx);
%Momentum equation
dvdt_c(j)=-v_nc(j)*dvdx-(1/gamma)*(dtdx+(t_nc(j)/rho_nc(j))*drhodx);
end
%Average
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 j=2:n-1
v_nc(j)=v_nc_old(j)+dvdt(j)*dt;
rho_nc(j)=rho_nc_old(j)+drhodt(j)*dt;
t_nc(j)=t_nc_old(j)+dtdt(j)*dt;
end
%Applying boundary condition
%Inlet
v_nc(1)=2*v_nc(2)-v_nc(3);
t_nc(1)=2*t_nc(2)-t_nc(3);
rho_nc(1)=2*rho_nc(2)-rho_nc(3);
%Outlet
v_nc(n)=2*v_nc(n-1)-v_nc(n-2);
t_nc(n)=2*t_nc(n-1)-t_nc(n-2);
rho_nc(n)=2*rho_nc(n-1)-rho_nc(n-2);
%calculating mass flow rate,mach number & pressure
mfl_nc=rho_nc.*v_nc.*a;
M_nc=v_nc./((t_nc).^0.5);
p_nc=rho_nc.*t_nc;
%Plotting mass flow rate at different time steps
figure(1)
if k==1
plot(x,mfl_nc,'*','color','b','markersize',10)
hold on
elseif k==50
plot(x,mfl_nc,'y','linewidth',1.5)
hold on
elseif k==100
plot(x,mfl_nc,'c','linewidth',1.5)
hold on
elseif k==200
plot(x,mfl_nc,'g','linewidth',1.5)
hold on
elseif k==600
plot(x,mfl_nc,'r','linewidth',1.5)
hold on
elseif k==800
plot(x,mfl_nc,'k','linewidth',1.5)
hold on
elseif k==1000
plot(x,mfl_nc,'--','color','bla','linewidth',1.5)
hold on
xlabel('Non_dimensional distance x')
ylabel('Mass flow rate')
legend('1 time step','50 time steps','100 time steps','200 time
steps',
'600 time steps','800 time steps','1000 time steps')
legend('Location','north')
title(['Non-Conversation analysis';'Mass flow rate at different time
steps'])
grid on
end
end
end

MAIN CODE:-In this code we have called the above two functions.And then the result is plotted
like density,velocity,temperature & velocity versus non-dimensional distance along
nozzle for both conservative form as well as for non-conservative form.

clear all
close all
clc
%Creating meshing & assiging Inputs
n=31;% number of grids points
x=linspace(0,3,n);
dx=x(2)-x(1);
gamma=1.4;
nt=1200; % number of time steps
%Calling the function for conservation form
tic
[rho_con,v_con,t_con,mfl_con,m_con,p_con]=Conservation(n,nt);
time_con=toc;
%Calling the function for Non_conservation form
tic
[rho_nc,v_nc,t_nc,mfl_nc,M_nc,p_nc]=Non_Conservation(n,nt);
time_nc=toc;
%Comparing mass flow rate for conservation and non conservation
figure(3)
plot(x,mfl_nc,'*','color','c','markersize',10)
hold on
plot(x,mfl_con,'color','r','linewidth',1)
xlabel('Non_dimensional distance x')
ylabel('Mass flow rate')
text1=sprintf('Conservation analysis time=%d',time_con);
text2=sprintf('Non_Conservation analysis time=%d',time_nc);
title=(['Comparsion of conservation & non conservation mass flow
rate';text1;text2])
legend('Conservation','Non_Conservation','location','northwest')
##%Plotting different parameters like density,velocity,temperature,mach
number & pressure
figure(4)
subplot(5,1,1)
plot(x,M_nc,'r','linewidth',1.5)
hold on
plot(x,m_con,'*','color','r','markersize',10)
ylabel('Mach number')
legend('Non-conservation','Conservation','location','northwest')
grid on
subplot(5,1,2)
plot(x,rho_nc,'r','linewidth',1.5)
hold on
plot(x,rho_con,'*','color','r','markersize',10)
ylabel('Density')
legend('Non-conservation','Conservation','location','northwest')
grid on
subplot(5,1,3)
plot(x,t_nc,'r','linewidth',1.5)
hold on
plot(x,t_con,'*','color','r','markersize',10)
ylabel('Temperature')
legend('Non-conservation','Conservation','location','northwest')
grid on
subplot(5,1,4)
plot(x,p_nc,'r','linewidth',1.5)
hold on
plot(x,p_con,'*','color','r','markersize',10)
ylabel('pressure')
legend('Non-conservation','Conservation','location','northwest')
grid on
subplot(5,1,5)
plot(x,v_nc,'r','linewidth',1.5)
hold on
plot(x,v_con,'*','color','r','markersize',10)
ylabel('Velocity')
legend('Non-conservation','Conservation','location','northwest')
title('Variation of Non-Dimensional flow parameters along x axis')
grid on

Results:-

1. Comparing the mass flow rate of different time steps between conservation &
non-conservation In order to figure out what should be the minimum number of cycles
for which the simulation should be run in order to get converge:-
Conservation:-

Non-Conservation:-

2. Comparing the normalized mass flow rate between the conservative forms and
non-conservative forms:-

Observation:-

1. In the above figure we can observe that non-conservation plot vary a lot like its
decreasing till throat area and then again increasing.This is because the control
volume is moving with the flow and hence there is a lot variation in the parameter.

3. Variation of Non-Dimensional flow parameters along x-axis:-

Observation:-
● In the above figure, we can observe that, plot are identical for conservation and
non-conservation, so we can say that mathematically conservation &
non-conservation is the same.