Challenge

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

SSugumar P6 Aug 2022

In this challenge you are going to simulate the isentropic flow through a quasi 1D subsonic-supersonic nozzle. You will derive both the conservation and non-conservation forms of the governing equations and sovle them using the MacCormack's technique. You need to determine the steady-state temperature distribution for the flow-field variables and investigate the difference between the two forms of governing equations by comparing their solutions.

You need to show the following plots inside your report
1. Steady-state distribution of primitive variables inside the nozzle
2. Time-wise variation of the primitive variables
3. Variation of Mass flow rate distribution inside the nozzle at different time steps during the time-marching process
4. Comparison of Normalized mass flow rate distributions of both forms

Let us consider a small control volume of the nozzle as shown below.

We assume that at a given section, where the cross section area is A, the flow properties are uniform across the section. Hence, although the area of the nozzle changes as a function of distance along the nozzle, X, and therefore in reality the flow field is two dimensional(the flow varies in the two dimensional x-y space), We make the assumption that the flow properties vary only with X; this is tantamount to asuming uniform flow properties across any given cross-section. Such flow is defined as Quasi one dimensional flow.

The governing continuty, momentum and energy equations in terms of non dimensional variables for this quassi one dimensional staedy and isentropc flow can be expressed respectively as:

Conservative Continuty equation:

∂(ρ′A′)∂t′+∂(ρ′A′V′)∂x′=0">∂(ρ'A')∂t'+∂(ρ'A'V')∂x'=0∂(ρ′A′)∂t′+∂(ρ′A′V′)∂x′=0

Conservative Momentum equation:

∂(ρ′A′V′)∂t′+∂[ρ′A′(V′)2+(1γ)ρ′A′]∂x′=(1γ)ρ′(∂A′∂x′)">∂(ρ'A'V')∂t'+∂[ρ'A'(V')2+(1γ)ρ'A']∂x'=(1γ)ρ'(∂A'∂x')∂(ρ′A′V′)∂t′+∂[ρ′A′(V′)2+(1γ)ρ′A′]∂x′=(1γ)ρ′(∂A′∂x′)

Conservative Energy equation:

∂[ρ′(e′γ-1+(γ2)(V′)2)A′]∂t′+∂[ρ′(e′γ-1+(γ2)(V′)2)V′A′+ρ′A′V′]∂x′=0">∂[ρ'(e'γ−1+(γ2)(V')2)A']∂t'+∂[ρ'(e'γ−1+(γ2)(V')2)V'A'+ρ'A'V']∂x'=0∂[ρ′(e′γ-1+(γ2)(V′)2)A′]∂t′+∂[ρ′(e′γ-1+(γ2)(V′)2)V′A′+ρ′A′V′]∂x′=0

For air at standard conditions, γ=1.4">γ=1.4γ=1.4

Non- Conservative Continuty equation:

∂ρ′∂t′=-ρ′(∂V′∂x′)-ρ′V′(∂(lnA′)∂x′)-V′(∂ρ′∂x′)">∂ρ'∂t'=−ρ'(∂V'∂x')−ρ'V'(∂(lnA')∂x')−V'(∂ρ'∂x')∂ρ′∂t′=-ρ′(∂V′∂x′)-ρ′V′(∂(lnA′)∂x′)-V′(∂ρ′∂x′)

Non- Conservative 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′))

Non- Conservative 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′]

These equations are solved by Macormack method using time marching process

Macormack Method:

Let f(x) is the equation to be solved.Then the following steps are followed

i. inital profile

Before the calculation, the initial profile for all the flow field variables which are velocity, density and temperature has to be determined. These are guess values but should be specified with more awareness on the problem.

ii. predector step

(∂fdt)p=f(i+1)-f(i)Δx">(∂fdt)p=f(i+1)−f(i)Δx(∂fdt)p=f(i+1)-f(i)Δx

iii. solution update

fp=fInitial+(∂f∂t)p⋅Δt">fp=fInitial+(∂f∂t)p⋅Δtfp=fInitial+(∂f∂t)p⋅Δt

fp=">fp=fp= Flow field variables calculated after the predictor step

fInitial=">fInitial=fInitial=flow field variable in the initial step

Δt=">Δt=Δt=time step size

iv. Corrector step

(∂fpdt)c=f(i)-f(i-1)Δx">(∂fpdt)c=f(i)−f(i−1)Δx(∂fpdt)c=f(i)-f(i-1)Δx

(∂fdt)avg=0.5((∂fdt)p+(∂fpdt)c)">(∂fdt)avg=0.5((∂fdt)p+(∂fpdt)c)(∂fdt)avg=0.5((∂fdt)p+(∂fpdt)c)

v. Final solution update

fsol=fInitial+(∂f∂t)avg⋅Δt">fsol=fInitial+(∂f∂t)avg⋅Δtfsol=fInitial+(∂f∂t)avg⋅Δt

`f_(sol) =`Final solution of the flow field variables for the first time step

Matlab Code with explaination:

close all
clear all
clc

L=3;                    %length of nozzle
n=[31 62 93];           %number of grid points
nt=1:1400;              %number of time steps
cfl = 0.5;              %CFL number
for i = 1:length(n)
x=linspace(0,L,n(i));
[mfr_nc,mfr_throat_nc]=non_conservative(L,n(i),nt,cfl);
[mfr_c,mfr_throat_c]=conservative(L,n(i),nt,cfl);
figure(7)
subplot(2,1,1)
plot(nt,mfr_throat_c,'b',nt,mfr_throat_nc,'r','linewidth',1)
xlabel('Time steps (nt)')
ylabel('Mass flow rate')
legend('Conservative form','Non-Conservative form')
title('Mass flow rate using "Non-Conservative form vs Conservative form" at throat point i=16')

subplot (2,1,2)
plot(x,mfr_c,'b',x,mfr_nc,'r','linewidth',1)
xlabel('Length of nozzle')
ylabel('Mass flow rate')
legend('Conservative form','Non-Conservative form')
title('Mass flow rate using "Non-Conservative form vs Conservative form" in the nozzle')

figure(8)
subplot (2,1,1)
plot(nt,mfr_throat_c,'linewidth',1)
hold on
xlabel('Time steps (nt)')
ylabel('Mass flow rate')
legend('@31 nodes','@62 nodes','@93 nodes')
title('Conservative form: Grid independency test at throat point i=16')
subplot (2,1,2)
plot(nt,mfr_throat_nc,'linewidth',1)
hold on
xlabel('Time steps (nt)')
ylabel('Mass flow rate')
legend('@31 nodes','@62 nodes','@93 nodes')
title('Non-Conservative form: Grid independency test at throat point i=16')
figure(9)
subplot (2,1,1)
plot(x,mfr_c,'linewidth',1)
hold on
xlabel('Length of nozzle')
ylabel('Mass flow rate')
legend('@31 nodes','@62 nodes','@93 nodes')
title('Conservative form: Grid independency test along the nozzle')
subplot (2,1,2)
plot(x,mfr_nc,'linewidth',1)
hold on
xlabel('Length of nozzle')
ylabel('Mass flow rate')
legend('@31 nodes','@62 nodes','@93 nodes')
title('Non-Conservative form: Grid independency test along the nozzle')

end

Conservative function:

function [mfr_c,mfr_throat_c]=conservative(L,n,nt,cfl)

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

% initial conditions rho, a, v, T, M, mfr, p
for ii=1:length(x)
    if (x(ii)<=0.5)
        rho(ii) = 1;
        T(ii) = 1;
    elseif (x(ii)>0.5) & (x(ii)<=1.5)
        rho(ii) = 1-0.366*(x(ii)-0.5);
        T(ii)= 1-0.167*(x(ii)-0.5);
    elseif (x(ii)>1.5) & (x(ii)<=3.5)
        rho(ii) = 0.634-0.3879*(x(ii)-1.5);
        T(ii)= 0.833-0.3507*(x(ii)-1.5);
    end
end
a= 1+2.2*(x-1.5).^2;     %area
g = 1.4;                 %gamma value
% solution vectors u1, u2, u3
u1 = rho.*a;
u2 = 0.59*ones(1,n);
v=u2./u1;
u3 = u1.*((T./(g-1))+(g/2).*v.^2);

M= v./(T.^0.5);          %mach number
mfr = rho.*a.*v;         %mass flow rate
p=rho.*T;                %pressure


for z=1:31               %to find dt value
u=(T(z))^(1/2);
dt(z) = cfl*dx/(u + v(z));
end
dt=min(dt);

for i = 1:length(nt)
    %storing initial values
    u1_old = u1;
    u2_old = u2;
    u3_old = u3;
   
    for j=2:n-1
   %flux vectors
        f1 = u2;
        f2 = ((u2.^2)./u1)+((g-1)/g)*(u3-((g/2)*(u2.^2)./u1));
        f3 = ((g.*u2.*u3)./u1) - ((g*(g-1)/2)*(u2.^3)./(u1.^2));
   %source term
        j1 = (1/g)*rho(j)*T(j)*((a(j+1)-a(j))/dx);        
   %predictor step
        du1_dt_p(j)= -(f1(j+1)-f1(j))/dx;           %continuity   
        du2_dt_p(j)= -(f2(j+1)-f2(j))/dx + j1;      %momentum
        du3_dt_p(j)= -(f3(j+1)-f3(j))/dx;           %energy
   %updating solution vectors          
        u1(j) = u1_old(j) + du1_dt_p(j)*dt;
        u2(j) = u2_old(j) + du2_dt_p(j)*dt;
        u3(j) = u3_old(j) + du3_dt_p(j)*dt;
   %updating primitives
        rho=u1./a;
        T=(g-1).*(u3./u1-(g/2).*(u2./u1).^2);  
   %updating flux terms
        f1=u2;
        f2=u2.^2./u1+((g-1)/g).*(u3-(g/2)*u2.^2./u1);
        f3=g.*u2.*u3./u1-0.5*g*(g-1).*u2.^3./u1.^2;
   %source term for corrector step
        j2 = (1/g)*(rho(j)*T(j)*(a(j)-a(j-1))/dx);    
   %corrector step
        du1_dt_c(j) = -(f1(j)-f1(j-1))/dx;            %continuity   
        du2_dt_c(j) = -(f2(j)-f2(j-1))/dx + j2;       %momentum
        du3_dt_c(j) = -(f3(j)-f3(j-1))/dx;            %energy
    end
    %average time derivatives
        du1_dt_avg = 0.5*(du1_dt_p + du1_dt_c);
        du2_dt_avg = 0.5*(du2_dt_p + du2_dt_c);
        du3_dt_avg = 0.5*(du3_dt_p + du3_dt_c);
    %final corected values
    for k = 2:n-1
        u1(k) = u1_old(k) + du1_dt_avg(k)*dt;
        u2(k) = u2_old(k) + du2_dt_avg(k)*dt;
        u3(k) = u3_old(k) + du3_dt_avg(k)*dt;
    end
    %entry boundary conditions
        u1(1) = a(1);
        u2(1) = 2*u2(2)-u2(3);
        u3(1) = u1(1)*(T(1)/(g-1) + (g*(u2(1)/u1(1))^2)/2);
    %exit boundary conditions
        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);

    %Final update of primitives
        rho = u1./a;
        v = u2./u1;
        T=(g-1).*(u3./u1-(g/2).*(u2./u1).^2);  
        p = rho.*T;
        M= v./T.^0.5;
        error(k) =  abs(v(k)-v(k-1));
        mfr= u2;
  %to find values at throat point
    T_throat(i) = T(16);
    rho_throat(i) = rho(16);
    v_throat(i) = v(16);
    M_throat(i) = M(16);
    p_throat(i) = p(16);
    error_throat(i)=error(16);
    mfr_throat(i)=mfr(16);  
    
if n==31  
    for ii=[1:50:201 700]
    if ii==i
       figure(4)
     plot(x,mfr,'linewidth',1)
     hold on
     title ('Conservative form - Variation of Mass flow rate distribution inside the nozzle at different time steps during the time-marching process')
     xlabel('Length of nozzle')
     ylabel('Mass flow rate')   
    elseif i==1400
     plot(x,mfr,'*r')
     legend('@1','@51','@101','@151','@201','@700','@1400') 
    end
end
end
  
end
if n==31
figure(5)
     plot(x,v,'b',x,p,'r',x,T,'k',x,rho,'m',x,M,'g',x,mfr,'c')
     ylim([-0.1 3.7])
     title (sprintf('Time step = %d',i))
     xlabel('Length of nozzle')
     ylabel('Non dimensional parameters')
     legend('Velocity','Pressure','Temperature','Density','Mach number','Mass flow rate','location','northwest')
     title('Conservative form - Steady-state distribution of primitive variables inside the nozzle')

figure (6)
plot(nt,rho_throat)  
hold on
plot(nt,v_throat,'linewidth',1)   
plot(nt,p_throat,'linewidth',1) 
plot(nt,T_throat,'linewidth',1)    
plot(nt,M_throat,'linewidth',1)
plot(nt,mfr_throat,'linewidth',1)
xlabel('Time steps (nt)')
ylabel('Non dimensional parameters')
legend('Density','Velocity','Pressure','Temperature','Mach number','Mass flow rate')
title('Conservative form - Time-wise variation of the primitive variables at throat point i=16')
end
mfr_c=mfr;
mfr_throat_c=mfr_throat;

end

Non-Conservative function:

function [mfr_nc,mfr_throat_nc]=non_conservative(L,n,nt,cfl)


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

% initial conditions rho, a, v, T, mfr, M, p
rho = 1-0.3146*x;
T= 1-0.2314*x;
v=(0.1+1.09*x).*T.^0.5;
a= 1+2.2*(x-1.5).^2;
M= v./(T.^0.5);
mfr = rho.*a.*v;
p=rho.*T;

g = 1.4;

for z=1:31
    u=(T(z))^(1/2);
    dt(z) = cfl*dx/(u + v(z));
end
dt=min(dt);

for i = 1:length(nt)
  %storing initial values
    rho_old = rho;
    T_old = T;
    v_old = v;
  for j=2:n-1
  %predictor step
    drho_dx = (rho(j+1)-rho(j))/dx;
    dv_dx = (v(j+1)-v(j))/dx;
    dT_dx = (T(j+1)-T(j))/dx;
    da_dx = (log(a(j+1))-log(a(j)))/dx;
    drho_dt_p(j) = -rho(j)*dv_dx - rho(j)*v(j)*da_dx - v(j)*drho_dx;        %continuty
    dv_dt_p(j) = -v(j)*dv_dx - (1/g)*(dT_dx + (T(j)/rho(j))*drho_dx);       %momentum
    dT_dt_p(j) = -v(j)*dT_dx - (g-1)*T(j)*(dv_dx + v(j)*da_dx);             %energy
  %updating primitives
    v(j) = v_old(j) + dv_dt_p(j)*dt;
    rho(j) = rho_old(j) + drho_dt_p(j)*dt;
    T(j) = T_old(j) + dT_dt_p(j)*dt;
  %corrector step
    drho_dx = (rho(j)-rho(j-1))/dx;
    dv_dx = (v(j)-v(j-1))/dx;
    dT_dx = (T(j)-T(j-1))/dx;
    da_dx = (log(a(j))-log(a(j-1)))/dx; 
    drho_dt_c(j) = -rho(j)*dv_dx - rho(j)*v(j)*da_dx - v(j)*drho_dx;        %continuity    
    dv_dt_c(j) = -v(j)*dv_dx - (1/g)*(dT_dx + (T(j)/rho(j))*drho_dx);       %momentum
    dT_dt_c(j) = -v(j)*dT_dx - (g-1)*T(j)*(dv_dx + v(j)*da_dx);             %energy
  end
  %average time derivatives
    drho_dt_avg = 0.5*(drho_dt_p + drho_dt_c);
    dv_dt_avg = 0.5*(dv_dt_p + dv_dt_c);
    dT_dt_avg = 0.5*(dT_dt_p + dT_dt_c);
  %final corected values
  for k = 2:n-1
    rho(k) = rho_old(k) + drho_dt_avg(k)*dt;
    v(k) = v_old(k) + dv_dt_avg(k)*dt;
    T(k) = T_old(k) + dT_dt_avg(k)*dt;
  end
  %entry boundary conditions
    v(1) = 2*v(2)-v(3);
  %exit boundary conditions
    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);
  %Final update of primitives
    p = rho.*T;
    M= v./T.^0.5;
    error(k) =  abs(v(k)-v(k-1));
    mfr= rho.*a.*v;
  %to find values at throat point
    T_throat(i) = T(16);
    rho_throat(i) = rho(16);
    v_throat(i) = v(16);
    M_throat(i) = M(16);
    p_throat(i) = p(16);
    error_throat(i)=error(16);
    mfr_throat(i)=mfr(16);
    if n==31
for ii=[1:50:201 700]
    if ii==i
       figure(1)
     plot(x,mfr,'linewidth',1)
     hold on
     title ('Non-Conservative form - Variation of Mass flow rate distribution inside the nozzle at different time steps during the time-marching process')
     xlabel('Length of nozzle')
     ylabel('Mass flow rate')   
    elseif i==1400
     plot(x,mfr,'*r')
     legend('@1','@51','@101','@151','@201','@700','@1400') 
    end
end
end
  
end
if n==31
figure(2)
     plot(x,v,'b',x,p,'r',x,T,'k',x,rho,'m',x,M,'g',x,mfr,'c')
     ylim([-0.1 3.7])
     title (sprintf('Time step = %d',i))
     xlabel('Length of nozzle')
     ylabel('Non dimensional parameters')
     legend('Velocity','Pressure','Temperature','Density','Mach number','Mass flow rate','location','northwest')
     title('Non-Conservative form - Steady-state distribution of primitive variables inside the nozzle')

figure (3)
plot(nt,rho_throat)  
hold on
plot(nt,v_throat,'linewidth',1)   
plot(nt,p_throat,'linewidth',1) 
plot(nt,T_throat,'linewidth',1)    
plot(nt,M_throat,'linewidth',1)
plot(nt,mfr_throat,'linewidth',1)
xlabel('Time steps (nt)')
ylabel('Non dimensional parameters')
legend('Density','Velocity','Pressure','Temperature','Mach number','Mass flow rate')
title('Non-Conservative form - Time-wise variation of the primitive variables at throat point i=16')
end
mfr_nc=mfr;
mfr_throat_nc=mfr_throat;
end

Explanation:

The purpose of the function is to calculate the values of
Mass flow rate, Pressure, Mach number, Density, velocity, Temperature, Density at throat, velocity at throat, Temperature at throat, Mass flow rate at throat, Pressure at throat, Mach number at throat by the known equations for all the grid points.
To determine the initial profile of Density & Temperature, for loop is used inputing the known conditions within an else-if command.
Initial values of area, throat, velocity & Pressure are calculated using formulas.
In order to simplify and avoid errors, equation is splitted into Solution Vector(U), Flux Vector(F) and Source Term(J).
Now Outer Time loop is computed for for loop to calculate all the grid points for various formulas inside the loop.
MacCarmack technique is used to solve the numerical solution of governing equations such as Continuity, Momentum & Energy Equation.
The technique implies two methods for solving namely Predictor method followed by Corrector method.
Predictor & Corrector method are assigned within the for loop individually for calculating values.
The calculated values are averaged in the next step and updating the solution with old values.
Boundary Conditions are applied within the for loop and subsequently used for calculating the Mass Flow Rate, Mach Number & Pressure and also calculating the values throat section.
Finally plotting is done for Distance(x) & Mass Flow Rate for differrent Timesteps.
Program is ended.