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

Objective : Numerical simulation of 1D supersonic nozzle flow using Macormack Method using conservative and non-conservative forms of equation

Assumptions:

1. Flow inside the nozzle is assumed to be isentropic.

2. The flow is considered Quasi-1D as properties vary along the x axis and not the y axis.

3. Compreesibe flow

4. No friction and heat transfer considerations.

5. Considerable pressure ration is maintained between choked pressure and back ressure to allow expansion of the flow and avoid shock waves.

Problem Description:

1. Area profile along the direction of the flow :

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

2. Initial non-dimnetionalized thermodynamic properties :

`rho = 1-0.3146*x`

`T = 1-0.2314*x`

`V = (0.1+1.09*x)^0.5`

Governing Equations:

1. Non-conservative form

• Continuity Equation : `(delrho)/(delt) = -rho*(delV)/(delx)-rho*V*(dellnA)/(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*(dellnA)/(delx))`

2. Conservative Form

• Continuity Equation : `(del(rhoA))/(delt) + (del(rhoVA))/(delx) = 0`
• Momentum Equation : `(del(rhoAV))/(delt)+(del(rhoAV^2+(1/gamma)PA))/(delx) = (1/gamma)*P*(delA)/(delx)`
• Energy Equation : `(del(rho*(T/(gamma-1)+(gamma/2)*V^2)*A))/(delt) + (del(rho*(T/(gamma-1)+(gamma/2)*V^2)*V*A+P*A*V))/(delt) = 0`

The equations were solved using Macormack Method.

Scripts:

Main script:

clear all
close all
clc

%Inputs
n = 31; %Number of nodes
x = linspace(0,3,n); %Mesh
dx = x(2) - x(1);
gamma = 1.4;

%Time steps
nt = 5000;
c1 = 0.5;
c2 = 0.5;
tol = 1e-6;
mass_tol1 = 1e-3;
mass_tol2 = 1e-2;

[rho1,v1,T1,total_time1,netmf1,err_v1,err_rho1,err_T1] = nonconserv(n,x,dx,gamma,nt,c1,tol,mass_tol1);

[rho2,v2,T2,total_time2,netmf2,err_v2,err_rho2,err_T2] = conserv(n,x,dx,gamma,nt,c2,tol,mass_tol2);

%Velocity plots
figure(1)
plot(x,v1,'b',x,v2,'r')
xlabel('Non-dimentional Distance')
ylabel('Non-dimentional Velocity')
legend('Non Conservative','Conservative')
grid on


%Density plots
figure(2)
plot(x,rho1,'b',x,rho2,'r')
xlabel('Non-dimentional Distance')
ylabel('Non-dimentional Density')
legend('Non Conservative','Conservative')
grid on

%Temperature plots
figure(3)
plot(x,T1,'b',x,T2,'r')
xlabel('Non-dimentional Distance')
ylabel('Non-dimentional Temperature')
legend('Non Conservative','Conservative')
grid on

%Net Mass Flow Rate Plot
figure(4)
subplot(2,1,1)
plot(netmf1)
xlabel('Time Step')
ylabel('Net Mass Flow Rate')
title('Non Conservative Form')
grid on
subplot(2,1,2)
plot(netmf2)
xlabel('Time Step')
ylabel('Net Mass Flow Rate')
title('Conservative Form')
grid on

%Error in Velocity with Time
figure(5)
subplot(2,1,1)
plot(err_v1)
xlabel('Time Step')
ylabel('Error in Velocity')
title('Non Conservative Form')
grid on
subplot(2,1,2)
plot(err_v2)
xlabel('Time Step')
ylabel('Error in Velocity')
title('Conservative Form')
grid on

%Error in density with Time
figure(6)
subplot(2,1,1)
plot(err_rho1)
axis([0 inf -0.01 0.02])
xlabel('Time Step')
ylabel('Error in Density')
title('Non Conservative Form')
grid on
subplot(2,1,2)
plot(err_rho2)
axis([0 inf -0.005 0.015])
xlabel('Time Step')
ylabel('Error in Density')
title('Conservative Form')
grid on

%Error in Temperature with Time
figure(7)
subplot(2,1,1)
plot(err_T1)
axis([0 inf -0.005 0.01])
xlabel('Time Step')
ylabel('Error in Temperature')
title('Non Conservative Form')
grid on
subplot(2,1,2)
plot(err_T2)
axis([0 inf -0.01 0.02])
xlabel('Time Step')
ylabel('Error in Temperature')
title('Conservative Form')
grid on

Conservative Form (function script):

function [rho,v,T,k,netmf,err_v,err_rho,err_T] = conserv(n,x,dx,gamma,nt,c,tol,mass_tol)

%Intial Profiles
rho = 1 - 0.3146*x; %Density
T = 1 - 0.2312*x; %Temperature
a = 1 + 2.2*(x-1.5).^2; %Area
v = 0.59./(rho.*a); %Velocity

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


% Outer Time Loop
for k = 1:nt
    
    % storing old values
    U1_old = U1;
    U2_old = U2;
    U3_old = U3;
    
    rho_old = U1_old./a;
    v_old = U2_old./(rho.*a);
    T_old = (gamma-1)*((U3_old./U1_old) - (gamma/2)*v.^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);
    
    for j = 2:n-1
        J2(j) = (1/gamma)*rho(j)*T(j)*(a(j+1) - a(j))/dx;
    end
    
    for i = 1:n
        deltat(i)=(c*(dx/((T(i)^0.5+v(i)))));
    end
    dt = min(deltat);
    
    % Predictor Method
    for j = 2:n-1
        %Continuity Equation
        dU1dt_p(j) = -(F1(j+1)-F1(j))/dx;
        
        %Momentum Equation
        dU2dt_p(j) = -(F2(j+1)-F2(j))/dx + J2(j);
        
        %Energy Equation
        dU3dt_p(j) = -(F3(j+1)-F3(j))/dx;
        
        %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;
    end
    
    %Calculating Primitive variables
    rho = U1./a;
    v = U2./(rho.*a);
    T = (gamma-1)*((U3./U1) - (gamma/2)*v.^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);
    
    for j = 2:n-1
        J2(j) = (1/gamma)*rho(j)*T(j)*(a(j) - a(j-1))/dx;
    end
    
    % Corrector Method
    for j = 2:n-1
        %Continuity Equation
        dU1dt_c(j) = -(F1(j)-F1(j-1))/dx;
        
        %Momentum Equation
        dU2dt_c(j) = -(F2(j)-F2(j-1))/dx + J2(j);
        
        %Energy Equation
        dU3dt_c(j) = -(F3(j)-F3(j-1))/dx;
    end
    
    %Computing average time derivative
    dU1dt = 0.5*(dU1dt_p + dU1dt_c);
    dU2dt = 0.5*(dU2dt_p + dU2dt_c);
    dU3dt = 0.5*(dU3dt_p + dU3dt_c);
    
    %Final Update
    for j = 2:n-1
        U1(j) = U1_old(j) + dU1dt(j)*dt;
        U2(j) = U2_old(j) + dU2dt(j)*dt;
        U3(j) = U3_old(j) + dU3dt(j)*dt;
    end
    
    %Applying Boundary Conditions
    %inlet
    U2(1) = 2*U2(2) - U2(3);
    
    %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);
    
    %Recalculating Primitive variables
    rho = U1./a;
    v = U2./(rho.*a);
    T = (gamma-1)*((U3./U1) - (gamma/2)*v.^2);
    
    %Mass Flow Rate
    inflow = rho(1)*v(1)*a(1);
    outflow = rho(n)*v(n)*a(n);
    netmf(k) = outflow-inflow;
    
    %Error Tolerance
    err_v(k) = max(v - v_old);
    err_rho(k) = max(rho - rho_old);
    err_T(k) = max(T - T_old);
    if(err_v(k) < tol && err_rho(k) < tol && err_T(k) < tol && netmf(k)<mass_tol)
        break;
    end
end
end

Non-Conservative Form (function script):

function [rho,v,T,k,netmf,err_v,err_rho,err_T] = nonconserv(n,x,dx,gamma,nt,c,tol, mass_tol)

%Intial Profiles
rho = (1 - 0.3146*x); %Density
T = (1 - 0.2314*x); %Temperature
v = (0.1 + (1.09*x)).*(T.^0.5); %Velocity

a = (1 + 2.2*(x-1.5).^2); %Area
% Outer Time Loop
for k = 1:nt
    
    dt = min(abs(c*dx./(realsqrt(T)+v)));
    
    %Copying old variables
    rho_old = rho;
    v_old = v;
    T_old = T;
    
    
    % Predictor Method
    for j = 2:n-1
        dvdx = (v(j+1) - v(j))/dx;
        drhodx = (rho(j+1) - rho(j))/dx;
        dlogadx = (log(a(j+1)) - log(a(j)))/dx;
        dTdx = (T(j+1) - T(j))/dx;
        
        %Continuity Equation
        drhodt_p(j) = -rho(j)*dvdx - rho(j)*v(j)*dlogadx - v(j)*drhodx;
        
        %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)*dlogadx);
        
        %Solution Update
        v(j) = v(j) + dvdt_p(j)*dt;
        rho(j) = rho(j) + drhodt_p(j)*dt;
        T(j) = T(j) + dTdt_p(j)*dt;
    end
    
    % Corrector Method
    for j = 2:n-1
        dvdx = (v(j) - v(j-1))/dx;
        drhodx = (rho(j) - rho(j-1))/dx;
        dlogadx = (log(a(j)) - log(a(j-1)))/dx;
        dTdx = (T(j) - T(j-1))/dx;
        
        %Continuity Equation
        drhodt_c(j) = -rho(j)*dvdx - rho(j)*v(j)*dlogadx - 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)*dlogadx);
    end
    
    %Computing average time derivative
    drhodt = 0.5*(drhodt_p + drhodt_c);
    dvdt = 0.5*(dvdt_p + dvdt_c);
    dTdt = 0.5*(dTdt_p + dTdt_c);
    
    %Final Update
    for j = 2:n-1
        v(j) = v_old(j) + dvdt(j)*dt;
        rho(j) = rho_old(j) + drhodt(j)*dt;
        T(j) = T_old(j) + dTdt(j)*dt;
    end
    
    %Applying 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);
    
    %Mass Flow Rate
    inflow = rho(1)*v(1)*a(1);
    outflow = rho(n)*v(n)*a(n);
    netmf(k) = (outflow-inflow);
    
    %Error Tolerance
    err_v(k) = max(v - v_old);
    err_rho(k) = max(rho - rho_old);
    err_T(k) = max(T - T_old);
    if(err_v(k) < tol && err_rho(k) < tol && err_T(k) < tol && netmf(k)<mass_tol)
        break;
    end
end
end

Result:

The simulation was run for a maximum of 5000 time-steps, with CFL=0.5 and n=31. The tolerance was kept to 1e-6, where the change in rho, v, and T had to fall under the tolerance. Net Mass Flow Rate tolerance has also been set at 1e-3 for non conservation form and 1e-2 for conservation form.

Non Conservation Form converged in 676 time steps, by converged it has reached a stable solution.

Conservation Form converged in 4515 time steps, by converged it has reached a stable solution.

Velocity Plots:

Density Plots:

Temperature Plot:

It has been observed that conservative form solves in less number of iterations than non-conservative form. However, the net mass flow rate hasn't stabilised in the current condition.

Net Mass Flow Rate :

Net Change in Velocity:

Net Change in Density:

Net Change in Temperature: