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

 The objective in this project is to write code to solve the 1D supersonic nozzle flow equations using the Macormack Method.

The governing equations for both conservative and nonconservative are solved and compared. Also, the iteration number and computational time are compared until the final profile is obtained. From the comparison, we can say that the conservative part is far better than the nonconservative form.

Here, we are using the MacCormack method for solving the 1D Super-sonic nozzle flow. Initially, the grid points are chosen as a minimum value and then the grid points are increased to see the changes in the solutions. At 1st, the grid is taken as 31 because this is the maximum value for the steady state to attain. At this points, the graphs are plotted and the result is given below. at the 2nd time, it is made at 51 which show the large variation.

The equations of the form are given below.

Non-conservative form

• Continuity Equation : 
• `∂ρ)/(∂t)=−ρ⋅(∂V)/(∂x)−ρ⋅V⋅(∂lnA)/(∂x)−V⋅(∂ρ)/(∂x)`
•  
•  
• Momentum Equation : ∂V∂t=-V⋅∂V∂x-(1γ)⋅(∂T∂x+(Tρ)⋅∂ρ∂x)∂V∂t=−V⋅∂V∂x−(1γ)⋅(∂T∂x+(Tρ)⋅∂ρ∂x)∂V∂t=-V⋅∂V∂x-(1γ)⋅(∂T∂x+(Tρ)⋅∂ρ∂x)
•  
•  
• Energy Equation : ∂T∂t=-V⋅∂T∂x-(γ-1)⋅T⋅(∂V∂x+V⋅∂lnA∂x)\">∂T∂t=−V⋅∂T∂x−(γ−1)⋅T⋅(∂V∂x+V⋅∂lnA∂x)∂T∂t=-V⋅∂T∂x-(γ-1)⋅T⋅(∂V∂x+V⋅∂lnA∂x)

2. Conservative Form

• Continuity Equation : ∂(ρA)∂t+∂(ρVA)∂x=0\">∂(ρA)∂t+∂(ρVA)∂x=0∂(ρA)∂t+∂(ρVA)∂x=0
•  
•  
• Momentum Equation : ∂(ρAV)∂t+∂(ρAV2+(1γ)PA)∂x=(1γ)⋅P⋅∂A∂x\">∂(ρAV)∂t+∂(ρAV2+(1γ)PA)∂x=(1γ)⋅P⋅∂A∂x∂(ρAV)∂t+∂(ρAV2+(1γ)PA)∂x=(1γ)⋅P⋅∂A∂x
•  
•  
• Energy Equation : ∂(ρ⋅(Tγ-1+(γ2)⋅V2)⋅A)∂t+∂(ρ⋅(Tγ-1+(γ2)⋅V2)⋅V⋅A+P⋅A⋅V)∂t=0\">∂(ρ⋅(Tγ−1+(γ2)⋅V2)⋅A)∂t+∂(ρ⋅(Tγ−1+(γ2)⋅V2)⋅V⋅A+P⋅A⋅V)∂t=0

----------------------------------------------------------------------------------

Nonconservative form

close all
clear all
clc

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

% Time step
dt = 0.01;
nt = 900;
gamma = 1.4;

% Initial condition
rho = 1 - 0.3146*x;
T = 1-.2314*x;    % T represents temperature
v = (0.1 + 1.09*x).*(T.^(.5));

% Area
a = 1 + 2.2*((x - 1.5).^2);

n_throat=find(a-1<0.001)

for k = 1:nt
    
    rho_old = rho;
    T_old = T;
    v_old = v;
    a_old=sqrt(T_old);
    massflow_old=rho_old.*v_old.*a;
  
    
    % Predictor step
    for j = 2:n-1
        
        drho_dx = (rho(j+1) - rho(j))/dx;
        dv_dx = (v(j+1) - v(j))/dx;
        dT_dx = (T(j+1) - T(j))/dx;
        dlog_a_dx = (log(a(j+1)) - log(a(j)))/dx;
        
        %contunity Equation
        drho_dt_p(j) = - rho(j)*dv_dx - rho(j)*v(j)*dlog_a_dx - v(j)*drho_dx;
       
        %Momentum equation
        dv_dt_p(j) = -v(j)*dv_dx - (dT_dx + (T(j)/rho(j))*drho_dx)/gamma;
       %energy Equation
        dT_dt_p(j) = -v(j)*dT_dx - (gamma-1)*T(j)*(dv_dx + v(j)*dlog_a_dx);
        
        
         % solution update 
        rho(j) = rho(j) + drho_dt_p(j)*dt;
        v(j) = v(j) + dv_dt_p(j)*dt;
        T(j) = T(j) + dT_dt_p(j)*dt;
    end
    
    % Corrector step
    for j = 2:n-1
        
        drho_dx = (rho(j) - rho(j-1))/dx;
        dv_dx = (v(j) - v(j-1))/dx;
        dT_dx = (T(j) - T(j-1))/dx;
        dlog_a_dx = (log(a(j)) - log(a(j-1)))/dx;
       
        %Contunity Equation
        drho_dt_c(j) = - rho(j)*dv_dx - rho(j)*v(j)*dlog_a_dx - v(j)*drho_dx;
      
        %Momentum Equation
        dv_dt_c(j) = -v(j)*dv_dx - (dT_dx + (T(j)/rho(j))*drho_dx)/gamma;
       
        %energy Equation
        dT_dt_c(j) = -v(j)*dT_dx - (gamma-1)*T(j)*(dv_dx + v(j)*dlog_a_dx);
              
    end
    %compute the average time derivatives
    drho_dt_avg = 0.5*(drho_dt_c + drho_dt_p);
    dv_dt_avg = 0.5*(dv_dt_c + dv_dt_p);
    dT_dt_avg = 0.5*(dT_dt_c + dT_dt_p);
    
    for i = 2:n-1
        %Final solution updates
        rho(i) = rho_old(i) + drho_dt_avg(i)*dt;
        v(i) = v_old(i) + dv_dt_avg(i)*dt;
        T(i) = T_old(i) + dT_dt_avg(i)*dt;
    end
    
    %Applying boundary condition
    
    % Inlet condition
    rho(1)=1;
    T(1)=1;
    v(1) = 2*v(2)-v(3);
    
    
    % Outlet condition
    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);
    
   
   
    %calculating
    %Mach number
    M=v./sqrt(T);
    %mass flow rate
    massflow=rho.*v.*a;
    p=rho.*T;
    
  %rho,T,P,M Throat Values
    rho_throat(k)=rho(n_throat);
    T_throat(k)=T(n_throat);
    p_throat(k)=p(n_throat);
    M_throat(k)=M(n_throat);
    
 
    %Check for time convergence 
    error=max(abs(massflow-massflow_old));

    
    TransientMassFlow=figure(1);
    hold on;
    
    if k== 1
        plot(x,massflow,\'--k\');
    else if k == 50
        plot(x,massflow,\'b\');
    else if k == 100
            plot(x,massflow,\'r\');
        else if k == 150
                plot(x,massflow,\'g\');
            else if k == 200
                    plot(x,massflow,\'b\');
                else if k == 700
                        plot(x,massflow,\'-.b\',\'linewidth\',1.25);
                        plot(x(end),massflow(end),\'k*\');
                    end 
                    end
                end
            end
        end
    end
    
end

title(\'Mass Flow Rates(Non-Conservative Form)\');
xlabel(\'x/L\');
ylabel(\'\\rhoVA/\\rho_{o}V_{o}A_{o}\');
legend(\'0\\Deltat\',\'50\\Deltat\',\'100\\Deltat\',\'150\\Deltat\',\'200\\Deltat\',\'700\\Deltat\');

grid on;

figure(2)
subplot(4,1,1);
plot(x,rho)
legend(\'Density ratio\');
ylabel(\'\\rho/\\rho_{o}\');
xlabel(\'X/L\');
grid on;

subplot(4,1,2)
plot(x,v,\'b\');
legend(\'velocity\');
ylabel(\'v/v{o}\');
xlabel(\'X/L\');
grid on;

subplot(4,1,3)
plot(x,T,\'b\')
legend(\'temperature\')
ylabel(\'T/T_{o}\');
xlabel(\'X/L\');
grid on;

figure(3)
plot(x,v,\'r\')
ylabel(\'non dimensionless velocity\')
xlabel(\'non dimensionless distance\')
legend(\'non conservative\')

figure(4)
plot(x,rho,\'g\')
ylabel(\'non dimensionless density\')
xlabel(\'non dimensionless distance\')
legend(\'non conservative\')
    

conservative 

close all
clear all
clc

% Inputs
n = 31;
x = linspace(0,3,n);
dx= x(2)-x(1);
rho = ones(1,n);
T = ones(1,n);
gamma = 1.4;
error=1e9;
error2=1e9;
tol=1e-5;

% Area
a = 1 + 2.2*((x - 1.5).^2);

n_throat=find(a-1<0.001)




for i = 1:n
    if x(i) <= 0.5
        rho(i) = 1;
        T(i) = 1;
    elseif x(i) > 0.5 && x(i)<= 1.5
        rho(i) = 1 - 0.366*(x(i) - 0.5);
        T(i) = 1 - 0.167*(x(i) -0.5);
    elseif x(i) > 1.5 && x(i) <= 3.5
        rho(i) = 0.634 - 0.3879*(x(i) - 1.5);
        T(i) = 0.833 - 0.3507*(x(i) -1.5);
    end
end


% velocity
U2 = .59; % For initial condition only
v = U2./(rho.*a);
p=rho(i).*T;


% Initial condition
U1 = rho.*a;
U2 = (rho.*a).*v;
U3 = rho.*((T./(gamma-1))+ (gamma/2)*v.^2).*a;

k2=0;


% Time marching 
nt = 1000;
dt = 0.0267;
massflow_old=rho.*v.*a;

for k = 1:nt
    
    U1_old = U1;
    U2_old = U2;
    U3_old = U3;
    
  
    
    
    
   %defining flux value
    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));

     % predictor step
    for i = 2:n-1
        
        J2 = (1/gamma)*rho(i)*T(i)*((a(i+1)-a(i))/dx);
        
        %contunity equation
        dU1_dt_p(i) = -(F1(i+1)-F1(i))/dx;
       %momentum equation
       dU2_dt_p(i) = -((F2(i+1)-F2(i))/dx) + J2;
       
       %energy equation
       dU3_dt_p(i) = -(F3(i+1)-F3(i))/dx;
        
       
        % solution update
        U1(i) = U1(i) + dU1_dt_p(i)*dt;
        U2(i) = U2(i) + dU2_dt_p(i)*dt;
        U3(i) = U3(i) + dU3_dt_p(i)*dt;
   %VALUES
         v=U2./U1;
        rho(i) = U1(i)/a(i);
        T(i) = (gamma-1)*((U3(i)/U1(i)) - ((gamma/2)*((U2(i)./U1(i)).^2)));
        p=rho(i);
    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 step
    for j = 2:n-1
        
        J2 = (1/gamma)*rho(j)*T(j)*((a(j)-a(j-1))/dx);
        
        dU1_dt_c(j) = -(F1(j)-F1(j-1))/dx;
        dU2_dt_c(j) = -((F2(j)-F2(j-1))/dx) + J2;
        dU3_dt_c(j) = -(F3(j)-F3(j-1))/dx;
        
    end
   
    
    % average solution update
    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);
    
    %final sol
    for r = 2:n-1
        U1(r) = U1_old(r) + dU1_dt(r)*dt;
        U2(r) = U2_old(r) + dU2_dt(r)*dt;
        U3(r) = U3_old(r) + dU3_dt(r)*dt;
    end
    
    rho = U1./a;
    v = U2./U1;
    T = (gamma-1)*(U3./U1 - (gamma/2)*(v.^2));
    
    % Boundary conditions
    
    % Inlet condition
    U1(1) = a(1);
    U2(1) = 2*U2(2)-U2(3);
    v(1) = U2(1)/U1(1);
    T(1) = 1;
    U3(1) = U1(1)*(T(1)/(gamma-1) + (gamma/2)*(v(1)^2));
    
    % Outlet condition
    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);
   
    rho(n) = U1(n)/a(n);
    v(n) = U2(n)/U1(n);
    T(n) = (gamma-1)*((U3(n)/U2(n))- (gamma/2)*(v(n)^2));
    p=rho(i).*T;
  
    %calculation of other variable
    M=v./sqrt(T);
      %mass flow rate
    massflow=rho.*v.*a;
    p=rho.*T;

   


    %Transient Throat Values
    rho_throatk(2)=rho(i).*(n_throat);
    T_throatk(2)=T(i).*(n_throat);
    p_throatk(2)=p(i).*(n_throat);
    M_throatk(2)=M(i).*(n_throat);
    
 
    %Check for time convergence 
    error2=max(abs(massflow-massflow_old));
    
     
    if k2 == 100
            plot(x,massflow_con,\'g\');
            else if k2 == 200
                    plot(x,massflow_con,\'r\');
                else if k2 == 700
                        plot(x,massflow_con,\'-.b\',\'linewidth\',1.25);
                end
            end
    end
end
plot(x,rho,\'r\')
title(\'Transient Non-dimensional Mass Flow Rates(Non-Conservative Form)\');
xlabel(\'x/L\');
ylabel(\'\\rhoVA/\\rho_{o}V_{o}A_{o}\');
legend(\'100\\Deltat\',\'200\\Deltat\',\'700\\Deltat\');
axis([0 3 0.5 0.8]);
grid minor;


subplot(4,1,2)
plot(x,v,\'b\');
legend(\'velocity\');
ylabel(\'v/v{o}\');
xlabel(\'X/L\');
grid on;

subplot(4,1,3)
plot(x,T,\'b\')
legend(\'temperature\')
ylabel(\'T/T_{o}\');
xlabel(\'X/L\');
grid on;

figure(2)
plot(x,v,\'r\')
ylabel(\'non dimensionless velocity\')
xlabel(\'non dimensionless distance\')
legend(\'conservative\')

figure(2)
subplot(4,1,1);
plot(x,rho)
legend(\'Density ratio\');
ylabel(\'\\rho/\\rho_{o}\');
xlabel(\'X/L\');
grid on;

figure(4)
plot(x,rho,\'g\')
ylabel(\'non dimensionless density\')
xlabel(\'non dimensionless distance\')
legend(\'conservative\')
    

Result:

The stability of the nonconservative form solver has been provided below. The flow rate is also high lighted at the inlet and an outlet portion. Also, we have found the all the possible values at the Mach number =1 at the throat.The simulation was run for a maximum of 900 time-steps, with n=31. The tolerance was kept to 1e-5, where the change in rho, v, and T had to fall under the tolerance. The error is considered pf 1e-9 for conservative form.

we can say, that conservative form solves in less number of iterations than non-conservative form. However, the net mass flow rate is in unstable conditions.

The transient Non-dimensional mass flow rates are given below;

The conservative form is given below;