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

The objective of this challenge is to perform simulation of a Quasi 1D Supersonic nozzle by Macormack method using conservative and non-conservative forms of governing equation

To find:

Density

Pressure

Temperature

Mach Number

Mass flow rate

Grid Dependency test for n = 31 and n = 61.

Number of Cycles at which the solution converges.

NOZZLE EXPLANATION:

The nozzle used here is the convergent-divergent nozzle.

The initial and boundary condition are provided by Neuman for boundary and Dirichlet for initial conditions.

The simulation is obtained for both the conservative form and nonconservative forms of governing equation.

CONSERVATIVE FORM VS NONCONSERVATIVE FORM:

The conservative and non-conservative forms od equations are mathematically the same. The difference comes when discretizing the equations. Notice that both the non-conservative and conservative forms come from the conservation equations (i.e. conservation of mass and conservation of momentum). So, you can have both a conservative form of the conservation equations and a non-conservative form of the conservation equations.

In the conservative form, the flux terms and force vectors are involved thus discretizing it is much heavy when compared to the neoconservative form.

NON-CONSERVATIVE FORM:

CONSERVATIVE FORM:

Continuity equation:

Momentum equation:

Energy equation:

Flux terms to Discretize the equation:

PRIMITIVE TERMS:

TIMESTEP CALCULATION:

MACORMACK METHOD:

To solve the set of governing equations,

1. Predictor method
2. Update the solution
3. Corrector method
4. Average of Predictor and corrector method.
5. Update the old values of the primitive terms.
6. Provide Boundary conditions.

FUNCTION CODE FOR NON-CONSERVATIVE FORM:

function [rho,v,t,Mass_Flow_Rate,Pressure_nc,M,x,rho_throat, v_throat, t_throat] = nonconservativenozzle(n)


x = linspace(0,3,n);
dx = x(2) - x(1);

rho = 1-(0.3146*x);
t = 1-(0.2314*x); % t = temperature
v = (0.1+1.09*x).*t.^0.5;
gamma = 1.4;
%area
a = 1+ 2.2*(x-1.5).^2;

 
C = 0.5; 
for j = 1:n
   
    dt(j) = C*(dx/(realsqrt(t(j))+v(j)));
end

dt = min(dt);

nt = 1400;

for z = 1:nt
    
rho_old = rho;
v_old = v;
t_old = t;
for j = 2:n-1
      
        dvdx = (v(j+1) - v(j))/dx;
        drhodx = (rho(j+1) - rho(j))/dx;
        dtdx = (t(j+1) - t(j))/dx;
        dlnadx = (log(a(j+1))-log(a(j)))/dx;
        
        
        
        %continuity equation
        A = -rho(j)*(dvdx); 
        B = -rho(j)*v(j)*(dlnadx);
        C = -v(j)*(drhodx);
        drhodt_p(j) = A+B+C;
        
         
        %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)*t(j)*(dvdx + v(j)*(dlnadx));
        
        
        %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
for j = 2:n-1
    
    
      
        dvdx = (v(j) - v(j-1))/dx;
        drhodx = (rho(j) - rho(j-1))/dx;
        dtdx = (t(j) -t(j-1))/dx;
        dlnadx = (log(a(j))-log(a(j-1)))/dx;
        
         %continuity equation
        drhodt_c(j) = -rho(j)*dvdx - rho(j)*v(j)*dlnadx - 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)*(t(j))*(dvdx + v(j)*(dlnadx));
        
        
        
     drhodt_avg(j) = 0.5*(drhodt_p(j) + drhodt_c(j));
     dvdt_avg(j) = 0.5*(dvdt_p(j) + dvdt_c(j));
     dtdt_avg(j) = 0.5*(dtdt_p(j) +dtdt_c(j));
     
        
end
for i = 2:n-1
     rho(i) = rho_old(i) + drhodt_avg(i)*dt;
     v(i) = v_old(i) + dvdt_avg(i)*dt;
     t(i) = t_old(i) + dtdt_avg(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);
     
     rho_throat(z) = rho(15);
     v_throat(z)  = v(15);
     t_throat(z) = t(15);
       
  
end

M = v./ ((t).^(0.5));
Mass_Flow_Rate = rho.*v.*a;
Pressure_nc = rho.*t;

end

FUNCTION CODE FOR CONSERVATIVE FORM:

function [rho1,V,T,Mass_Flow_rate,Pressure,Mach_Number,x,rho1_throat, V_throat, T_throat] = conservativenozzle(n)


x = linspace(0,3,n);
dx = x(2) - x(1);
A = 1+2.2*(x-1.5).^2;
gamma = 1.4;
%initial Conditions
for k = 1:n
    
     if x(k)>= 0 && x(k)<= 0.5
        rho1(k) = 1;
        T(k) = 1;
       elseif x(k)>0.5 && x(k)<= 1.5
        rho1(k) =  1 - 0.366*(x(k) - 0.5);
        T(k) = 1 - 0.167*(x(k) - 0.5);
        elseif x(k)> 1.5 && x(k)<= 3.5
        rho1(k) =  0.634 - 0.3879*(x(k) - 1.5);
        T(k) = 0.833 - 0.3507*(x(k) - 1.5);
     end
end

V = 0.59./(rho1.*A);

U1 = rho1.*A;
U2 = rho1.*A.*V;
U3 = rho1.*((T/(gamma-1))+((gamma)/2)*V.^2).*A;





nt = 1400;
for z = 1:nt
    for j = 1:n
   
    dt(j) = 0.5*(dx/(realsqrt(T(j))+V(j)));
    end
 dt = min(dt);
    
    
    U1_old = U1;
    U2_old = U2;
    U3_old = U3;

    F1_p = U2;
    F2_p = (U2.^2./U1) + ((gamma - 1)/gamma)*(U3 -(gamma/2)*(U2.^2./U1));
    F3_p = (gamma.*U2.*U3./U1) - ((gamma*(gamma-1)/2)*U2.^3./U1.^2);
 for j = 2:n-1
    
         
         
        dU1dt_p(j) = -(F1_p(j+1) - F1_p(j))/dx;
        J2_p(j) = (1/gamma)*rho1(j).*T(j).*(A(j+1) -A(j))/dx;
        dU2dt_p(j) = -(F2_p(j+1) - F2_p(j))/dx + J2_p(j);
        dU3dt_p(j) = -(F3_p(j+1) - F3_p(j))/dx;
        
        
        U1(j) = U1(j) + dU1dt_p(j)*dt;
        U2(j) = U2(j) + dU2dt_p(j)*dt;
        U3(j) = U3(j) + dU3dt_p(j)*dt;
        
 end
 
 rho1 = U1./A;
 
 T = (gamma-1)*((U3./U1)-((gamma/2).*(U2./U1).^2));
 
          
F1_c = U2;
F2_c = (U2.^2./U1) + ((gamma - 1)/gamma)*(U3 -(gamma/2)*(U2.^2./U1));
F3_c = (gamma.*U2.*U3./U1) - ((gamma*(gamma-1)/2)*U2.^3./U1.^2);

     
for j = 2:n-1
               
        dU1dt_c(j) = -(F1_c(j) - F1_c(j-1))/dx;
        J2_c(j) = (1/gamma).*rho1(j).*T(j).*(A(j) -A(j-1))/dx;
        dU2dt_c(j) = -(F2_c(j) - F2_c(j-1))/dx + J2_c(j);
        dU3dt_c(j) = -(F3_c(j) - F3_c(j-1))/dx;
        
        
        
        
end
         dU1dt_avg = 0.5*(dU1dt_c + dU1dt_p);
        dU2dt_avg = 0.5*(dU2dt_c + dU2dt_p);
        dU3dt_avg = 0.5*(dU3dt_c + dU3dt_p);
        
 for j = 2:n-1
        U1(j) = U1_old(j) + dU1dt_avg(j)*dt;
        U2(j) = U2_old(j) + dU2dt_avg(j)*dt;
        U3(j) = U3_old(j) + dU3dt_avg(j)*dt;
 end
   
  
     
    
    
      %BOUNDARY CONDITIONS FOR INLET
 
    U2(1) = 2*U2(2) - U2(3);
    V = U2./U1;
   U3(1) = rho1(1)*((T(1)/(gamma-1))+((gamma)/2)*V(1)^2).*A(1);
    
    
    %BOUNDARY CONDITIONS FOR 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);

    
    
    %PRIMITIVE TERMS
     rho1 = U1./A;
     V = U2./U1;
     T = (gamma-1)*((U3./U1)-(gamma/2).*(U2./U1).^2);
     
        rho1_throat(z) = rho1(15);
     V_throat(z)  = V(15);
     T_throat(z) = T(15);
end

Mach_Number = V./ ((T).^(0.5));
Mass_Flow_rate = rho1.*V.*A;
Pressure = rho1.*T;
end

MAIN CODE FOR COMPARING THE PRIMITIVE TERMS:

clear all
close all
clc

[rho1,V,T,Mass_Flow_rate,Pressure,Mach_Number,x]= conservativenozzle(31);

[rho,v,t,Mass_Flow_Rate,Pressure_nc,M,x] = nonconservativenozzle(31);



figure(1)
plot(x,v,x,V,'*')
xlabel('Length of Nozzle')
ylabel('Velocity')
legend('Non-Conservative form' , 'Conservative form')
figure(2)
plot(x,rho,x,rho1,'*')
xlabel('Length of Nozzle')
ylabel('Density')
legend('Non-Conservative form' , 'Conservative form')
figure(3)
plot(x,t,x,T,'*')
xlabel('Length of Nozzle')
ylabel('Temperature')
legend('Non-Conservative form' , 'Conservative form')
figure(4)
plot(x,M,x,Mach_Number,'*')
xlabel('Length of Nozzle')
ylabel('Mach Number')
legend('Non-Conservative form' , 'Conservative form')
figure(5)
plot(x,Pressure_nc,x,Pressure,'*')
xlabel('Length of Nozzle')
ylabel('Pressure')
legend('Non-Conservative form' , 'Conservative form')
figure(6)
plot(x,Mass_Flow_Rate,x,Mass_Flow_rate)
hold on
%exact solution rho'v'a' = 0.5790*ones(1,31);
y = 0.5790*ones(1,31);
plot(x,y, '--k')
xlabel('Length of Nozzle')
ylabel('Mass flow rate')
legend('Non-Conservative form' , 'Conservative form','Exact Solution')

MAIN CODE FOR TIME STEP VARIATION ANALYSIS AT THROAT:

clear all
close all
clc


[rho,v,t,Mass_Flow_Rate,Pressure_nc,M,x,rho_throat,v_throat,t_throat] = nonconservativenozzle(31);
[rho1,V,T,Mass_Flow_rate,Pressure,Mach_Number,x, rho1_throat,V_throat,T_throat]= conservativenozzle(31);


figure(1)
subplot(3,1,1)
plot([1:1400], rho_throat)
xlabel('No of tiem steps')
ylabel('rho/rho_0')
title('Timewise variation of Density, Temp & Velocity at throat(i=15) for non-conservative form')
subplot(3,1,2)
plot([1:1400], v_throat)
xlabel('No of tiem steps')
ylabel('v/v_0')
subplot(3,1,3)
plot([1:1400], t_throat)
xlabel('No of tiem steps')
ylabel('t/t_0')



figure(2)

subplot(3,1,1)
plot([1:1400], rho1_throat)
xlabel('No of tiem steps')
ylabel('rho/rho_0')
title('Timewise variation of Density, Temp & Velocity at throat(i=15) for Conservative form')
subplot(3,1,2)
plot([1:1400], V_throat)
xlabel('No of tiem steps')
ylabel('V/V_0')
subplot(3,1,3)
plot([1:1400], T_throat)
xlabel('No of tiem steps')
ylabel('T/T_0')

MAIN CODE FOR GRID DEPENDENCY TEST ON BOTH FORMS:

clear all
close all
clc



[rho,v,t,Mass_Flow_Rate,Pressure_nc,M,x] = nonconservativenozzle(31);
[rho_c,v1,t1,Mass_Flow_Rate1,Pressure_nc1,M1,x1] = nonconservativenozzle(61);


[rho1,V,T,Mass_Flow_rate,Pressure,Mach_Number,x1_1]= conservativenozzle(31);
[rho1_2,V1,T1,Mass_Flow_rate1,Pressure1,Mach_Number1,x1_2]= conservativenozzle(61);



figure(1)
plot(x,v)
hold on
plot(x1,v1,'*')
plot(x,t)
hold on
plot(x1,t1,'*')
hold on
plot(x,rho)
hold on
plot(x1,rho_c,'*')
grid on
title('Grid dependency test on non-conservative form')
legend('Grid points n = 31', 'Grid points n = 61')
xlabel('Length of Nozzle')
ylabel('Density, Temp , Velocity')
grid on



figure(2)
plot(x1_1,rho1)
hold on
plot(x1_2,rho1_2,'*')
hold on
plot(x1_1,T)
hold on
plot(x1_2,T1,'*')
hold on
plot(x1_1,V)
hold on
plot(x1_2,V1,'*')
grid on
title('Grid dependency test on Conservative form')
legend('Grid points n = 31', 'Grid points n = 61')
xlabel('Length of Nozzle')
ylabel('Density, Temp , Velocity')
grid on

RESULTS FOR PRIMITIVE TERMS:

Density ratio:

Density decreases across the area of the nozzle.

Velocity ratio:

Velocity increases across the area of the nozzle.

Temperature ratio:

Temperature decreases across the area of the nozzle.

Pressure ratio:

Pressure decreases across the area of the nozzle.

Mach Number:

Mach number is greater than 1 in the outlet and less than 1 in the inlet, the subsonic-supersonic nozzle.

MASS FLOW RATE:

The mass in conservative form is more conserved when compared to the non-conservative form.

RESULTS FOR TIME STEP VARIATION ANALYSIS AT THROAT:

For Conservative Form:

For Non-conservative form:

We could notice that after 400-500 time steps the flow attains a steady-state for both conservative and neoconservative form.

RESULTS FOR GRID DEPENDENCY TEST ON BOTH FORMS:

For Conservative Form:

The Density, temperature(r) and Velocity(b) are represented in a single plot.

For Non-conservative form:

The Density, temperature(b) and Velocity(r) are represented in a single plot

Hereby we can say that the solution produced is grid-independent. We obtain the same values even after doubling the grid points.

CONCLUSIONS:

• The conservation form yields a better mass flow distribution. The conservation
  form simply does a better job of conserving mass.
• The nonconservation form leads to smaller residuals. The amount by which the
  residuals decay is often used as an index of "quality" of the numerical algorithm.
  In this sense, the nonconservation form does a better job.
• There is no clear superiority of either form in terms of the accuracy of the results.
  The nonconservation form seems to produce slightly more accurate results for
  the primitive variables and the conservation form seems to produce slightly
  more accurate