BlogFor BusinessTalk to a counsellor
CoursesIIT ProgramsPlacementsLog in
BlogFor Business
Talk to a counsellor
Project

CFD || Simulation of a 1D Super-sonic Nozzle Flow using Macormack Method

JJayandran S Rao3 Aug 2019

OBJECTIVE

Write code solve the 1D supersonic nozzle flow equations using the Macormack Method.

  • Implement both the conservative and non-conservative forms of the governing equations
  • Perform grid dependence test

OVERVIEW

Conside steady,isentropic flow through a convergent-divergent nozzle as shown in Fig.The flow expands isentropically to supersonic speeds at the nozzle exit. The flow is locally subsonic in the convergent section of the nozzle, sonic at the throat and supersonic at the divergent section.We assume thatthe flow properties vary only with x; ie uniform flow properties any given c/s.Such flow is called quasi-1D flow.

  • Maccromack's Technique

Maccromack,s method is an explicit finite-difference technique which is second order in both space and time.Here a two-step predictor corrector sequence is usedwith forward differences on the predictor and with reward difference on the corrector.It is much easier to apply because there is no need to evaluate the second time derivatives.(For more info check "Computational Fluid Dynamics The Basic with application" by John D Anderson)

  • Governing Equations
  1. Conservative Equations

`(partial(rho'A'))/(partialt')+(partial(rho'A'V'))/(partialx') = 0` (Continuity Eqn)

`(partial(rho'A'V'))/(partialt')+(partial(rho'A'V'^2+(p'A')/gamma))/(partialx')=(p')/gamma(partialA')/(partialx')` (Momentum Eqn)

`(partial[rho'((e')/(gamma-1)+gamma/2V'^2)A'])/(partialt')+(partial[rho'((e')/(gamma-1)+gamma/2V'^2)A'V'+p'A'V'])/(partialx')=0` (Energy Eqn)

2. Non-conservative Equations

`(partialrho')/(partialt')=-rho'(partialV')/(partialx')-rho'V'partial/(partialx')(lnA')-V'(partialrho')/(partialx')` (Continuity Eqn)`

`(partialV')/(partialt')=-V'(partialV')/(partialx')-1/gamma((partialT)/(partialx')+(T')/(rho')(partialrho')/(partialx'))` (Momentum Eqn)`

`(partialT)/(partialt')=-V'(partialT')/(partialx')-(gamma-1)T'[(partialV')/(partialx')+V'partial/(partialx')(lnA')]` (Energy Eqn)`

Here the flow field are non-dimensionalized, where the flow variable are referenced to their reservior values.

`T'=T/T_0 ; rho'=rho/rho_0 ; x'=x/L ; V'=V/a_0 ; t' = t/(L/a_0) ; A' = A/A^*`

PROGRAM

  • Main function
%%%%%%%%%%%MAIN FUNCTION%%%%%%%%%%%

clear all
close all
clc

n=31;
x=linspace(0,3,n);
dx=x(2)-x(1);
gama=1.4;
a=1+2.2*(x-1.5).^2;
nt=1000;
C = 0.5; % Courant Number
throat = (n-1)/2;

%Solution for NON_CONSERVATIVE FROM
tic
[m,M,v,t,rho,p,t_throat,rho_throat,p_throat,M_throat]=Non_conservative(n,x,dx,gama,a,nt,C,throat);
timencons=toc;

%Plotting for non conservative
figure(2)
subplot(4,1,1);
plot(1:nt,rho_throat,\'b\');
xlabel(\'No.of Time steps\');
ylabel(\'rho/rho_{0}\')
legend(\'Density ratio\')
grid on;
grid minor;

subplot(4,1,2);
plot(1:nt,t_throat,\'r\');
xlabel(\'No.of Time steps\');
ylabel(\'T/T_{0}\');
legend(\'Temperature ratio\')
grid on;
grid minor;

subplot(4,1,3);
plot(1:nt,p_throat,\'m\');
xlabel(\'No.of Time steps\');
ylabel(\'P/P_{0}\');
legend(\'Pressur ratio\')
grid on;
grid minor;

subplot(4,1,4);
plot(1:nt,M_throat,\'g\');
xlabel(\'No.of Time steps\');
ylabel(\'M\');
legend(\'Mach number\')
grid on;
grid minor;



%Solution for conservative form
tic
[m_c,M_c,v_c,t_c,rho_c,p_c,t_throat_c,rho_throat_c,p_throat_c,M_throat_c]=conservative(n,x,dx,gama,a,nt,C,throat);
timecons=toc;

%Plotting for conservative
figure(4)
subplot(4,1,1);
plot(1:nt,rho_throat_c,\'b\');
xlabel(\'No.of Time steps\');
ylabel(\'rho/rho_{0}\')
legend(\'Density ratio\')
grid on;
grid minor;

subplot(4,1,2);
plot(1:nt,t_throat_c,\'r\');
xlabel(\'No.of Time steps\');
ylabel(\'T/T_{0}\');
legend(\'Temperature ratio\')
grid on;
grid minor;

subplot(4,1,3);
plot(1:nt,p_throat_c,\'m\');
xlabel(\'No.of Time steps\');
ylabel(\'P/P_{0}\');
legend(\'Pressur ratio\')
grid on;
grid minor;

subplot(4,1,4);
plot(1:nt,M_throat_c,\'g\');
xlabel(\'No.of Time steps\');
ylabel(\'M\');
legend(\'Mach number\')
grid on;
grid minor;

%Comparison of mass flow rate plot

figure(5)
plot(x,m,\'r\')
hold on
plot(x,m_c,\'g\')
hold on
line([0 3], [0.579 0.579])
xlabel(\'Non dimensional distance through nozzle(x/L)\')
ylabel(\'Non dimensional mass flow\')
title(\'Comparision of mass flow rate\')
legend(\'Non Conservative Soln\',\'Conservative Soln\',\'Exact Soln\')
  • Conservative equation
function [m_c,M_c,v_c,t_c,rho_c,p_c,t_throat_c,rho_throat_c,p_throat_c,M_throat_c]=conservative(n,x,dx,gama,a,nt,C,throat)

%Intial profile
for i=1:length(x)
    if x(i) >=0 && x(i)< 0.5
        rho_c(i) = 1;
        t_c(i) = 1;
    
    else if x(i) >=0.5 && x(i)< 1.5
            rho_c(i) = 1.0 - 0.366*(x(i)-0.5);
            t_c(i) = 1.0 - 0.167*(x(i)-0.5);
       
        else if x(i) >=1.5 && x(i)<= 3.5
                rho_c(i) = 0.634 - 0.3879*(x(i)-1.5);
                t_c(i) = 0.833 - 0.3507*(x(i)-1.5);
            end
        end
    end
end
v_c = 0.59./(rho_c.*a);
p_c = rho_c.*t_c;

%initializing solution vector
u1 = rho_c.*a;
u2 = rho_c.*a.*v_c;
u3 = rho_c.*a.*((t_c./(gama-1)) + ((gama/2)*(v_c.^2)));

% Time Iteration Loop
for k = 1:nt
    
    % Saving the old solution vectors:
    u1_old = u1;
    u2_old = u2;
    u3_old = u3;
    
    %Time-step control:
    dt = min(C.*dx./(sqrt(t_c) + v_c));
    
    % Defining the flux vectors:
    f1 = u2;
    f2 = (u2.^2./u1) + ((gama-1)/gama)*(u3 - (gama*u2.^2)./(2*u1));
    f3 = ((gama*u2.*u3)./u1) - (gama*(gama-1)/2)*(u2.^3./u1.^2);
    
for i=2:n-1
        
        %Source Term
        j2(i)=(1/gama).*rho_c(i).*t_c(i).*((a(i+1)-a(i))/dx);
        
        %predictor step
        du1dt_p(i) = -(f1(i+1) - f1(i))/dx;
        du2dt_p(i) = -(f2(i+1) - f2(i))/dx + j2(i);
        du3dt_p(i) = -(f3(i+1) - f3(i))/dx;
        
        %Update
        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
%Primitrive value and flux vector update
        rho_c = u1./a;
        v_c = u2./u1;
        t_c = (gama -1)*((u3./u1) - ((gama/2)*v_c.^2));
        p_c = rho_c.*t_c;
        
        f1 = u2;
        f2 = (u2.^2./u1) + ((gama-1)/gama)*(u3 - (gama*u2.^2)./(2*u1));
        f3 = ((gama*u2.*u3)./u1) - (gama*(gama-1)/2)*(u2.^3./u1.^2);
        
        %Corrector step
        
        for i=2:n-1
        
        %Source Term
         j2(i) = (1/gama)*rho_c(i).*t_c(i)*((a(i)-a(i-1))/dx);
        
        
        du1dt_c(i) = -(f1(i) - f1(i-1))/dx;
        du2dt_c(i) = -(f2(i) - f2(i-1))/dx + j2(i);
        du3dt_c(i) = -(f3(i) - f3(i-1))/dx;
        
        end
        
        %Average value
        
        du1dt=0.5*(du1dt_p+du1dt_c);
        du2dt=0.5*(du2dt_p+du2dt_c);
        du3dt=0.5*(du3dt_p+du3dt_c);
        
        %FInal solution
        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
        
        %Boundary condition
        %u1(1)=rho_c(1)*a(1);
        u2(1)= 2*u2(2) - u2(3);
        %u3(1)=u1(1)*(t_c(1)/(gama-1) + (gama*(v_c.^2))/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);
         
        %Final solution of primitive variables
        rho_c = u1./a;
        v_c = u2./u1;
        t_c = (gama -1).*((u3./u1) - ((gama/2).*v_c.^2));
        p_c = rho_c.*t_c;
        m_c= rho_c.*a.*v_c;
        M_c=v_c./sqrt(t_c);
        
        %Capturing throat values
        
        t_throat_c(k)=t_c(throat);
        p_throat_c(k)=p_c(throat);
        rho_throat_c(k)=rho_c(throat);
        M_throat_c(k)=M_c(throat);
        
        %MASS FLOW RATE AT DIFFERENT TIME OF THROAT
    
    figure(3)
    hold on;
    if k==100
        plot(x,m_c,\'b\')
    else if k==200
            plot(x,m_c,\'r\')
        else if k==700
                plot(x,m_c,\'g\')            
            end
        end
    end
end
    title(\'Non dimensional mass flow at different time(CONSERVATIVE)\');
    xlabel(\'Non dimensional distance through nozzle(x/L)\');
    ylabel(\'Non dimensional mass flow\');
    legend(\'100 dt\',\'200 dt\',\'700 dt\');
    grid on;
    grid minor;
end
  • Non-conservative equation
function [m,M,v,t,rho,p,t_throat,rho_throat,p_throat,M_throat]=Non_conservative(n,x,dx,gama,a,nt,C,throat)

%intial Profile
rho=1-0.3146*x;
t=1-0.2314*x;
v=(0.1+1.09*x).*t.^0.5;

%MACROMEK METHOD

for k=1:nt
    
    %SAving intial value
    rho_old=rho;
    t_old=t;
    v_old=v;
    
    %Time control
    dt=min((C.*dx)./(sqrt(t)+v)); 

    
    for i=2:n-1
    
        dvdx=(v(i+1)-v(i))/dx;
        drhodx=(rho(i+1)-rho(i))/dx;
        dlogadx=(log(a(i+1))-log(a(i)))/dx;
        dtdx=(t(i+1)-t(i))/dx;

        %Predictor Method
        drhodt_p(i)=-rho(i)*dvdx -rho(i)*v(i)*dlogadx -v(i)*drhodx;
        dvdt_p(i)=-v(i)*dvdx -((1/gama)*(dtdx + (t(i)/rho(i))*drhodx));
        dtdt_p(i)=-v(i)*dtdx -((gama-1)*t(i)*(dvdx + v(i)*dlogadx));

        %Update
        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 Method
    for i=2:n-1
        
        
        dvdx=(v(i)-v(i-1))/dx;
        drhodx=(rho(i)-rho(i-1))/dx;
        dlogadx=(log(a(i))-log(a(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/gama)*(dtdx+(t(i)/rho(i))*drhodx));
        dtdt_c(i)=-v(i)*dtdx - ((gama-1)*t(i)*(dvdx+v(i)*dlogadx));
    end
    %Average value
    drhodt=0.5*(drhodt_p + drhodt_c);
    dvdt=0.5*(dvdt_p + dvdt_c);
    dtdt=0.5*(dtdt_p + dtdt_c);
   
    %final Solution 
    for j=2:n-1
    
        
        rho(j)=rho_old(j)+drhodt(j)*dt;
        v(j)=v_old(j)+dvdt(j)*dt;
        t(j)=t_old(j)+dtdt(j)*dt;
    end
    
    %Apply BC
    
    %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);
    
    p=rho.*t;
    m=rho.*a.*v;
    M=v./sqrt(t);
    
    %Capturing throat values
    
    t_throat(k)=t(throat);
    p_throat(k)=p(throat);
    rho_throat(k)=rho(throat);
    M_throat(k)=M(throat);
    
    %MASS FLOW RATE AT DIFFERENT TIME OF THROAT
    
    figure(1)
    hold on;
    if k==50
        plot(x,m,\'b\')
    else if k==100
            plot(x,m,\'r\')
        else if k==150
                plot(x,m,\'g\')
            else if k==200
                    plot(x,m,\'m\')
                else if k==700
                        plot(x,m,\'-\')
                    end
                end
            end
        end
    end
end
    title(\'Non dimensional mass flow at different time(NON CONSERVATIVE)\');
    xlabel(\'Non dimensional distance through nozzle(x/L)\');
    ylabel(\'Non dimensional mass flow\');
    legend(\'50 dt\',\'100 dt\',\'150 dt\',\'200 dt\',\'700 dt\');
    grid on;
    grid minor;
end

RESULT

  • Non-conservative

The above fig represent the time-wise variation of density,temperature,pressure and mach number at the throat

  • Conservative

The above fig represent the time-wise variation of density,temperature,pressure and mach number at the throat

Time taken for non-conservative=5.48 sec

Time taken for conservative=6.67 sec

  • Grid independency test

The above values are at the throat of the nozzle

CONCLUSION

  1. The conservative form does a better job in mass flow distribution.
  2. The conservative form requires more time compared to non-conservative form;this is due to decoding the primitive variables from flux variable which is not required while using non conservative form.
  3. The number of grid points n=31 is enough to produce valid accurate result and when increased further doesn't produce a large variation in the result compared to exact results.
Your portfolio, next

Ship work like this
in your own course.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.