Simulation of a 1D Super-sonic nozzle flow simulation using MacCormack Method

AIM: Simulation of a 1D Supersonic nozzle flow using MacCormack Method for conservative form and non-conservative form.

THEORY: The objective of project is to solve 1D Super-sonic nozzle flow using numerical method.The numerical method which is used is MacCormack Method which is an explicit finite difference technique which provides second order accuracy solution.

For simulation of 1D supersonic nozzle flow, we will be solving two forms of equation using the MacCormack method:

• Conservation form
• Non-conservation form.

Assumption - Flow property varies only in x-direction and keeping y-direction as a constant.

CONSERVATION FORM: Conservation form refers to an arrangement of an equation or system of equations, usually representing a hyperbolic, that emphasizes that a property represented is conserved.

Governing equation for Conservation form:

Continuity equation:

`(delU_1)/(dt')=-(delF_1)/(dx')`

Momentum equation:

`(delU_2)/(dt')=-(delF_2)/(dx')+J_2`

Energy equation

`(delU_3)/(dt')=-(delF_3)/(dx')`

Where,

Flux term:

`F_1=rho'A'V'`

`F_2=rho'A'V'^2+i/gammap'A'`

`F_3=rho'(t/(gamma-1)+gamma/2V'^2)V'A'+rho'A'V'`

Solution vector:

`U_1=rho'A'`

`U_2=rho'A'V'`

`U_3=rho'(t/(gamma-1)+gamma/2V'^2)A'`

Source term:-

`J_2=1/gammap(delA')/(delx')`

Mesh Formation: In this steps the mesh is generated by forming grids along x-direction of nozzle as shown

Initial profile:-In this step we assign initial profile for different parameters

For `0<=x'<=0.5 {(rho'=1.0), (T'=1.0):}`

For `0.5<=x'<=1.5 {(rho'=1.0-0.366(x'-0.5)), (T'=1.0-0.167(x'-0.5)):}`

For `1.5<=x'<=3.5 {(rho'=0.634-0.3879(x'-1.5)), (T'=0.833-0.3507(x'-1.5)):}`

`nu'=U_2=(0.59)/(rho'A'^2) ;p=rho'T';M=(nu')/(sqrtT')`

Discretization: The process of converting PDE equation into algebraic equation is basically discretization. 

• Step for the Mccormack method to solve it.

MacCormack Method

• Predicator Step: In this step we discretize the above governing equation by applying first order forward differencing scheme to R.HS.
• Solution Update: In this step we update the solution which got from predictorstep
• Calculating Primitive Variable: The solution which we got from the above step is for source term,solution vector & flux term so,in this step we calculate useful parameters.
• Corrector Step: In this step we discretize the above governing equation by applying first order backward differencing scheme to R.HS.
• Average Time Step: In this step we average the solution which we got from the predicator step & corrector step. Final solution update:- In this,step we get required solution

Boundary conditions:

Inlet

`u_1(1)=ρ(1)*a(1)`

`u_2(1)=2*u_2(2)-u_2(3)`

`u_3(1)=u_1(1)*[(t(1))/(γ-1)+(gamma/2)*v(1)^2]`

`v(1)=u_2(1)/u1(1)`

Outlet

`u_1(n)=2*u_1(n-1)-u_1(n-2)`

`u_2(n)=2*u_2(n-1)-u_2(n-2)`

`nbsp;`u_3(n)=2*u_3(n-1)-u_3(n-2)`

FUNCTION CODE CONSERVATION :

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.

Governing equation for Conservation form:

Continuity equation:-

`(delrho')/(delt')=-rho'(delV')/(delx')-rho'V'(del(lnA)')/(delx')-V'(delrho')/(delx')`

Momentum equation:

`(delV')/(delt')=-V'(delV')/(delx')-1/gamma((delT')/(delx')-(T')/(rho')(delrho')/(delx'))`

Energy equation:- 

`(delT')/(delt')=-V'(delT')/(delx')-(gamma-1)T'[(delV')/(delx')+V'(del(lnA'))/(delx')]`

Mesh Formation: In this steps the mesh is generated by forming grids along x-direction of nozzle as shown

Initial profile: In this step we assign initial profile for different parameter

`ρ=1-0.3146x`

`v=(0.1+1.09x)t^0.5`

`a=1+2.2(x-1.5)^2`

`t=1-0.2314x`

Discretization: The process of converting PDE equation into algebraic equation is basically discretization. 

• Step for the Mccormack method to solve it.

MacCormack Method

• Predicator Step: In this step we discretize the above governing equation by applying first order forward differencing scheme to R.HS.
• Solution Update: In this step we update the solution which got from predictorstep
• Calculating Primitive Variable: The solution which we got from the above step is for source term,solution vector & flux term so,in this step we calculate useful parameters.
• Corrector Step: In this step we discretize the above governing equation by applying first order backward differencing scheme to R.HS.
• Average Time Step: In this step we average the solution which we got from the predicator step & corrector step. Final solution update:- In this,step we get required solution

Boundary conditions: 

Inlet:

`v(1)=2*v(2)-v(3)`

`t(1)=2t(2)-t(3)`

`rho(1)=2rho(2)-rho(3)`

Outlet:

`v(n)=2v(n-1)-v(n-2)`

`t(n)=2t(n-1)-t(n-2)`

`ρ(n)=2*ρ(n-1)-ρ(n-2)`

FUNCTION CODE NON CONSERVATION :

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

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_con,'*','color','c','markersize',10)
hold on
plot(x,mfl_nc,'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:

Observations: Observing the above plots we can conclude that the minimum number of cycles required for convergence is 700 cycles because after that results are the same as we can see for 1000 cycles.

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

Observation:

• We can observe that,non-conservation plot varies a lot like its decreasing till throat area & then again increasing. This is because the control volume is moving with the flow and hence there is a lot of variation in the parameter.
• While conservation plot does not vary that much like non-conservation plot, this is because in this case control volume is fixed and that's why there is not much variation in the parameter.
• Conservation analysis time is 18.3204 sec 
• Non-conservation analysis time is 15.004 sec,so we can say that non-conservation analysis method is faster then the conservation analysis ,this is because in conservation analysis there are some extra calculations that need to be done to calculate some primitive terms, so that it takes extra time.

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

Observation:

We can observe that plots are identical for conservation and non-conservation,so we can say that mathematically conservation & nonconservation are the same.