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

In this assignment, unsteady flow through the 1-D convergent-divergent nozzle is explored using MacCormak's technique.

The transient flow is solved using finite difference methods on the non-dimensionalized governing equation's conservative and non-conservative forms.

Separate functions are created for conservative and non-conservative methods, and they are called in the primary function to evaluate respective quantities. 

close all
clear all
clc

L=3;                      % Length of the nozzle
n=31;                     % no of grid points
gamma=1.4;                % gamma of the fluid
C=0.5;                    % Courant number
nt=1400;                  % number of time steps

[rho_SNC,T_SNC,V_SNC,rho_throatNC,T_throatNC,V_throatNC]=nonconservative_flow(L,n,gamma,C,nt);
[rho_SC,T_SC,V_SC,rho_throatC,T_throatC,V_throatC]=conservative_flow(L,n,gamma,C,nt);

mass_throatNC=(rho_throatNC.*V_throatNC.*T_throatNC);
mass_throatC=(rho_throatC.*V_throatC.*T_throatC);

x=linspace(0,L,n); 
plot(1:nt,mass_throatNC,1:nt,mass_throatC);
xlabel(' distance along the nozzle ');
ylabel(' dimensionless massflow rate ');
title(' transient throat massflow rate distribution for both methods ');
legend(' Non-conservative ', ' Conservative ');

 The function for evaluating the properties using the non-conservative method is as follows.

function [rho_SNC,T_SNC,V_SNC,rho_throatNC,T_throatNC,V_throatNC]=nonconservative_flow(L,n,gamma,C,nt)

% rho_SC is the steady-state density after nt time steps
% T_SC is the steady-state temperature after nt time steps
% V_SC is the steady-state velocity after nt time steps

%n is no of grid points
x=linspace(0,L,n);        % domain discretization
dx=x(2)-x(1);             % grid length

%C is Courant number
%nt is the number of time steps

% Initial conditions for dimensionless quantities using nonconservative
rho=1-0.3146*x;           % density
T=1-0.2314*x;             % temperature
V=(0.1+1.09*x).*sqrt(T);  % velocity
A=1+2.2*(x-1.5).^2;       % area


for j=1:nt                % time marching
    
    rho_old=rho;
    V_old=V;
    T_old=T;
    
    for i=2:n-1
        
        deltat(i-1)=C*(dx/(sqrt(T(i))+V(i)));
        
    end
    
    dt=min(deltat);    % finding the minimum timestep dt
    
    % Predictor step
    for i=2:n-1
        
        dVdx=(V(i+1)-V(i))/dx;
        dlogAdx=(log(A(i+1))-log(A(i)))/dx;
        drhodx=(rho(i+1)-rho(i))/dx;
        dTdx=(T(i+1)-T(i))/dx;
        
        drhodt_p(i)=-rho(i)*dVdx-rho(i)*V(i)*dlogAdx-V(i)*drhodx;
        dVdt_p(i)=-V(i)*dVdx-(1/gamma)*(dTdx+(T(i)/rho(i))*drhodx);
        dTdt_p(i)=-V(i)*dTdx-(gamma-1)*T(i)*(dVdx+V(i)*dlogAdx);
            
    end
    
    for i=2:n-1
        
        rho(i)=rho(i)+drhodt_p(i)*dt;
        V(i)=V(i)+dVdt_p(i)*dt;
        T(i)=T(i)+dTdt_p(i)*dt;
        
    end
    
    % Corrector step
    for i=2:n-1
        
        dVdx=(V(i)-V(i-1))/dx;
        dlogAdx=(log(A(i))-log(A(i-1)))/dx;
        drhodx=(rho(i)-rho(i-1))/dx;
        dTdx=(T(i)-T(i-1))/dx;
        
        drhodt_c(i)=-rho(i)*dVdx-rho(i)*V(i)*dlogAdx-V(i)*drhodx;
        dVdt_c(i)=-V(i)*dVdx-(1/gamma)*(dTdx+(T(i)/rho(i))*drhodx);
        dTdt_c(i)=-V(i)*dTdx-(gamma-1)*T(i)*(dVdx+V(i)*dlogAdx);
            
    end
    
    for i=2:n-1 
        
        % evaluating the properties at the next time step
        
        rho(i)=rho_old(i)+0.5*(drhodt_p(i)+drhodt_c(i))*dt;
        V(i)=V_old(i)+0.5*(dVdt_p(i)+dVdt_c(i))*dt;
        T(i)=T_old(i)+0.5*(dTdt_p(i)+dTdt_c(i))*dt;
        
    end
    
    % Boundary conditions for the first and last node
    V(1)=2*V(2)-V(3);
    rho(1)=1;
    T(1)=1;
    
    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); 
    
    % transient throat properties
    rho_throatNC(j)=rho((n+1)/2); 
    T_throatNC(j)=T((n+1)/2);
    V_throatNC(j)=V((n+1)/2);
    
    rho_SNC=rho;
    V_SNC=V;
    T_SNC=T;
    
end

end

 The function for evaluating the properties using the conservative method is as follows.

function [rho_SC,T_SC,V_SC,rho_throatC,T_throatC,V_throatC]=conservative_flow(L,n,gamma,C,nt)

% rho_SNC is the steady-state density after nt time steps
% T_SNC is the steady-state temperature after nt time steps
% V_SNC is the steady-state velocity after nt time steps

%n is no of grid points
x=linspace(0,L,n);        % domain discretization
dx=x(2)-x(1);             % grid length

%C is Courant number
%nt is the number of time steps

% Initial conditions and dimensionless quantities
A=1+2.2*(x-1.5).^2;       % area

for i=1:n
    % Initial conditions
    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)
        rho(i)=0.634-0.3879*(x(i)-1.5);
        T(i)=0.833-0.3507*(x(i)-1.5);
    end
end

% Initial values
for i=1:n
        V(i)=0.59/(rho(i)*A(i));
        U1(i)=rho(i)*A(i);
        U2(i)=rho(i)*A(i)*V(i);
        U3(i)=rho(i)*A(i)*((T(i)/(gamma-1))+(gamma/2)*V(i)^2);
end

for j=1:nt
    
    U1_old=U1;
    U2_old=U2;
    U3_old=U3;
    
    for i=2:n-1
        deltat(i-1)=C*(dx/(sqrt(T(i))+V(i))); 
    end
    
    dt=min(deltat);
    
    for i=1:n
        F1(i)=U2(i);
        F2(i)=(U2(i)^2/U1(i))+(gamma-1)*(U3(i)-(gamma/2)*U2(i)^2/U1(i))/gamma;
        F3(i)=gamma*(U2(i)*U3(i)/U1(i))-((gamma*(gamma-1)*U2(i)^3)/(2*U1(i)^2));   
    end
    
    % Predictor step
    for i=2:n-1
        dF1dx=(F1(i+1)-F1(i))/dx;
        dF2dx=(F2(i+1)-F2(i))/dx;
        dF3dx=(F3(i+1)-F3(i))/dx;
        J2=rho(i)*T(i)*(A(i+1)-A(i))/(gamma*dx);
        
        dU1dt_p(i)=-dF1dx;
        dU2dt_p(i)=-dF2dx+J2;
        dU3dt_p(i)=-dF3dx;      
    end
    
    for i=2:n-1
        U1(i)=U1(i)+dU1dt_p(i)*dt;
        U2(i)=U2(i)+dU2dt_p(i)*dt;
        U3(i)=U3(i)+dU3dt_p(i)*dt;
    end
    
    for i=1:n
        rho(i)=U1(i)/A(i);
        V(i)=U2(i)/U1(i);
        T(i)=(gamma-1)*((U3(i)/U1(i))-(gamma/2)*V(i)^2);
    end
    
    for i=1:n
        F1(i)=U2(i);
        F2(i)=(U2(i)^2/U1(i))+(gamma-1)*(U3(i)-(gamma/2)*U2(i)^2/U1(i))/gamma;
        F3(i)=gamma*(U2(i)*U3(i)/U1(i))-((gamma*(gamma-1)*U2(i)^3)/(2*U1(i)^2));   
    end
    
    % Corrector step
    for i=2:n-1
        dF1dx=(F1(i)-F1(i-1))/dx;
        dF2dx=(F2(i)-F2(i-1))/dx;
        dF3dx=(F3(i)-F3(i-1))/dx;
        J2=rho(i)*T(i)*(A(i)-A(i-1))/(gamma*dx);
        
        dU1dt_c(i)=-dF1dx;
        dU2dt_c(i)=-dF2dx+J2;
        dU3dt_c(i)=-dF3dx;      
    end
    
    for i=2:n-1
        U1(i)=U1_old(i)+0.5*(dU1dt_p(i)+dU1dt_c(i))*dt;
        U2(i)=U2_old(i)+0.5*(dU2dt_p(i)+dU2dt_c(i))*dt;
        U3(i)=U3_old(i)+0.5*(dU3dt_p(i)+dU3dt_c(i))*dt;
    end
    
    U1(1)=rho(1)*A(1);
    U2(1)=2*U2(2)-U2(3);
    U3(1)=U1(1)*((T(1)/(gamma-1))+(gamma/2)*V(1)^2);
    
    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);
    
    for i=1:n
        rho(i)=U1(i)/A(i);
        V(i)=U2(i)/U1(i);
        T(i)=(gamma-1)*((U3(i)/U1(i))-(gamma/2)*V(i)^2);
    end 
    
    % Transient throat values
    rho_throatC(j)=rho((n+1)/2);
    T_throatC(j)=T((n+1)/2);
    V_throatC(j)=V((n+1)/2);
    
    rho_SC=rho;
    T_SC=T;
    V_SC=V;
       
end

end

The corresponding plots for n=31 obtained are as follows.

Both methods evaluate the steady-state results correctly to a good approximation. The non-conservative method shows more oscillations than the conservative method in determining the steady-state values. 

Grid independence test results are as follows.