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

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

AAbdul Lateef2 Feb 2020
MATLAB

OBJECTIVE: Simulation of a 1D super-sonic nozzle flow simulation using the Macormack method

we implement both the conservative and non-conservative forms of the governing equations

and write code for solving 1D supersonic nozzle flow equation using the Macormack method.

we also perform the grid-dependent test.

MACORMACK METHOD:

Macormack method is used to provide second order accuracy without involvement of the second derivative, which helps in reducing the complexity.

Mach number: It is the ratio of the speed of the body to the sound speed.

If mach number < 1 = subsonic flow,

If mach number > 1 = supersonic flow.

Governing equations for 1D nonconservative flow in the non-dimensional form:

Continuity:

Momentum:

Energy:

In order to solve these above equations in matlab we need to follow the following steps:

  1. Create a mesh to solve the equation,
  2. Discretization,
  3. Solve the discretized equations using the boundary conditions.

we are using the Macormack method to solve the equations:

Algorithm to solve

Time marching problem

Initial conditions :

MATLAB CODE:

function nonconservative:

function[rho,t,v,P,mach_num,rhont,Pnt,tnt,mach_numnt,massraten1dt,massraten50dt,massraten100dt,massraten200dt,massraten500dt,massraten700dt,massraten1400dt] = fun_nonconser(L,n,x,dx,gamma,k)

    %initial profile
    %density
    rho = 1-0.3146*x;
    %t = temp
    t = 1-0.2314*x;
    % velocity
    v = (0.1+1.09*x).*t.^0.5;
    % Area
    a = 1+2.2*(x-1.5).^2;
    %time steps

    nt = 1400;
    dt = 0.01;
    throat = find(a==1);
    massraten1dt = zeros(n);
    massraten200dt = zeros(n);
    massraten400dt = zeros(n);
    massraten600dt = zeros(n);
    massraten800dt = zeros(n);
    massraten1000dt = zeros(n);
    massraten1200dt = zeros(n);
    massraten1400dt = zeros(n);
    
    %time loop

for k=1:nt
    
    
    rho_old = rho;
    v_old = v;
    t_old = t;
    
    %predictior terms
for j=2:n-1
           
    dvdx = (v(j+1)-v(j))/dx;
    drhodx = (rho(j+1)-rho(j))/dx;
    dtdx = (t(j+1)-t(j))/dx;
    dlogadx = (log(a(j+1))-log(a(j)))/dx;
    
    %countinuty equation predictor step
    
    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 
   
    rho(j) = rho(j) + drhodt_p(j)*dt;
    v(j) = v(j) + dvdt_p(j)*dt;
    t(j) = t(j) + dtdT_p(j)*dt;
    
end
    
     %corrector terms
for j=2:n-1
           
    dvdx = (v(j)-v(j-1))/dx;
    drhodx = (rho(j)-rho(j-1))/dx;
    dtdx = (t(j)-t(j-1))/dx;
    dlogadx = (log(a(j))-log(a(j-1)))/dx;
    
    %countinuty equation corrector step
    
    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 the avg of time derivative
    
    
    drhodt = 0.5*(drhodt_p + drhodt_c);
    
    dvdt = 0.5*(dvdt_p + dvdt_c);
    
    dtdT = 0.5*(dtdT_p + dtdT_c);
    
    %final solution update
for i=2:n-1
    rho(i) = rho_old(i) + drhodt(i)*dt; 
    v(i) = v_old(i) + dvdt(i)*dt;  
    t(i) = t_old(i) + dtdT(i)*dt;  
   
end
    
   %apply boundary conditions
   %inlet
    v(1) = 2*v(2)-v(3);

   %outlet
   rho(n) = 2*rho(n-1)-rho(n-2);
   v(n) = 2*v(n-1)-v(n-2);
   t(n) = 2*t(n-1)-t(n-2);
  
   %pressure value 
   P = rho.*t;
   mach_num = v./sqrt(t);
   massflow = rho.*v.*a;
   %throat values
   mach_numnt(k)= mach_num(throat);
   Pnt(k) = P(throat);
   rhont(k) = rho(throat);
   tnt(k) = t(throat);
   


  %saving mass floe rate at different time steps

 if k == 1 
    massraten1dt = massflow ;   
    elseif k == 50
    massraten50dt = massflow ;
    elseif k == 100
      massraten100dt = massflow ; 
    elseif k == 200
      massraten200dt = massflow ;  
    elseif k == 500
      massraten500dt = massflow ;
    elseif k == 700
      massraten700dt = massflow ;
    elseif k == 1400
      massraten1400dt = massflow ;
       
 end

figure(2)

   plot(x,v)
   xlabel(\'steps\')
   ylabel(\'velocity ratio\')
   title(sprintf(\'variation of non conser velocity for every itteration %d\',k ));
   drawnow
   F(k) = getframe(gcf);

end

    video = VideoWriter(\'onednozzlenonconder\',\'Uncompressed AVI\');
    open(video)
    writeVideo(video,F);
    close(video)
end

In order to see the minimum number of cycles required for convergence :

we require around 800 to 900 to get the convergence and the smooth curve as shown.

Governing equations for 1D conservative flow in the non-dimensional form:

Continuity:

Momentum:

Energy:

In conservative we use several constant terms to simplify the calculation

Let us define the elements of the solutions vector U, the flux vector F, and the source term J as follows.

we get these equations

To obtain the primitive variables (p, V, T,p, etc.), we must decode the elements U1, U2, and U3 as follows. From the definitions of U1 , U2 , and U3 given above, we have

PURE FORM OF THE FLUX TERMS. Consider the flux terms :

we get,

The flow properties a! the downstream, supersonic outflow boundary are obtained by linear extrapolation from the two adjacent internal points. If N denotes the gnd pomt at the outflow boundary, then

Boundary conditions:

In this also similar steps follow as above using macormak method

MATLAB CODE:

FUNCTION CONSERVATIVE:

function[rhoc,tc,vc,Pc,mach_numc,rhoct,Pct,tct,mach_numct,massrate1dt,massrate50dt,massrate100dt,massrate200dt,massrate500dt,massrate700dt,massrate1400dt] = fun_conser(L,n,x,dx,gamma,k)

%initial profile
%density variation 
for i = 1:n
    
   if (x(i)>=0 && 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

    % Area
    a = 1+2.2*(x-1.5).^2;
    % velocity
    v = 0.59./(rho.*a);
    %pressure 
    p = rho.*t;
    %time steps

    c = 0.5;

    dt = min(c.*(dx./((t.^0.5)+v)));
    
    nt = 1400;

    %dt = min(dt_grid);
    throat = find(a==1);

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

    %time loop

for  k=1:nt
    
    
    U1_old = U1;
    U2_old = U2;
    U3_old = U3;
   
    %defining flux 
    F1 = U2;  
    term = (((gamma-1)/gamma)*(U3-((gamma/2)*((U2.^2)./U1))));
    F2 = ((U2.^2)./U1)+ term;
    F3 = ((gamma*U2.*U3)./U1)-(gamma*(gamma-1)*(U2.^3))./(2*U1.^2);
    
    
     %predictior terms
    for j=2:n-1
       
    dF1dx(j) = (F1(j+1)-F1(j))/dx;
    dF2dx(j) = (F2(j+1)-F2(j))/dx;
    dF3dx(j) = (F3(j+1)-F3(j))/dx;
    dadx(j) = (a(j+1)-a(j))/dx;
    J2(j) = (1/gamma)*(rho(j)*t(j))*dadx(j);

    
    %countinuty equation predictor step
    
    dU1dt_p(j) = -dF1dx(j);
    
    %momentum equation
    
    dU2dt_p(j) = -dF2dx(j) + J2(j);
    
    %energy equation
    
    dU3dT_p(j) = -dF3dx(j);
    
   %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
    
   %updating rho v t values after predictor
    
    rhoc = U1./a;
    vc = U2./U1;
    tc = (gamma-1)*((U3./U1)-((gamma*vc.^2)/2));
    pc = rhoc.*tc;
    
    %updating flux
    
    F1 = U2;  
    F2 = ((U2.^2)./U1)+(((gamma-1)/gamma)*(U3-(gamma*U2.^2)./(2*U1)));
    F3 = ((gamma*U2.*U3)./U1)-(gamma*(gamma-1)*(U2.^3))./(2*U1.^2);
    
    
     %corrector terms
    for j=2:n-1
           
    dF1dx(j) = (F1(j)-F1(j-1))/dx;
    dF2dx(j) = (F2(j)-F2(j-1))/dx;
    dF3dx(j) = (F3(j)-F3(j-1))/dx;
    dadx(j) = (a(j)-a(j-1))/dx;
    J2(j) = (1/gamma)*(rhoc(j)*tc(j))*dadx(j);
    
    %countinuty equation corrector step
    
    dU1dt_c(j) = -dF1dx(j);
    
    %momentum equation
    
    dU2dt_c(j) = -dF2dx(j) + J2(j);
    
    %energy equation
    
    dU3dT_c(j) = -dF3dx(j);
    
 
    end
   
    %computing the avg of time derivative
    
    
    dU1dt = 0.5*(dU1dt_p + dU1dt_c);
    
    dU2dt = 0.5*(dU2dt_p + dU2dt_c);
    
    dU3dT = 0.5*(dU3dT_p + dU3dT_c);
    
    
    
    %final solution update
    for i=2:n-1
        
    U1(i) = U1_old(i) + dU1dt(i)*dt; 
    U2(i) = U2_old(i) + dU2dt(i)*dt;  
    U3(i) = U3_old(i) + dU3dT(i)*dt;  
   
    end
  
    %apply boundary conditions
    %inlet
    U(1) = rho(1).*a(1);
    U2(1) = 2*U2(2)-U2(3);
    U(3) = U1(1)*((t(1)/(gamma-1))+((gamma/2)*v(1)^2));
    
  
    %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);
    
   
    
    rhoc = U1./a;
    tc = (gamma-1)*(U3./U1-(gamma/2)*(U2./U1).^2);
    vc = U2./U1;
    mach_numc = vc./sqrt(tc);
    Pc = rhoc.*tc;
    massflow = rhoc.*vc.*a;
    %throat values
     mach_numct(k)= mach_numc(throat);
     Pct(k) = Pc(throat);
     rhoct(k) = rhoc(throat);
     tct(k) = tc(throat);
     
     
     if k == 1 
    massrate1dt = massflow ;   
    elseif k == 50
    massrate50dt = massflow ;
    elseif k == 100
      massrate100dt = massflow ; 
    elseif k == 200
      massrate200dt = massflow ;  
    elseif k == 500
      massrate500dt = massflow ;
    elseif k == 700
      massrate700dt = massflow ;
    elseif k == 1400
      massrate1400dt = massflow ; 
    
    end
  
   figure(1)
   plot(x,vc)
   xlabel(\'steps\')
   ylabel(\'velocity ratio\')
   title(sprintf(\'variation of conser velocity for every itteration %d\',k ));
   drawnow
    F(k) = getframe(gcf);
     
     
end
  
video = VideoWriter(\'onednozzleconser\',\'Uncompressed AVI\');
open(video)
writeVideo(video,F);
close(video)

end

In order to see the minimum number of cycles required for convergence :

It may take around 800 to 900 cycles to get the convergence and smooth curve as shown.

MAIN CODE:

clear all
close all
clc
%1d conservative
%inputs
L = 3;
n = 61;
x = linspace(0,3,n);
dx = x(2)-x(1);
gamma = 1.4;

k = 1400;
tic;
 [rho,t,v,P,mach_num,rhont,Pnt,tnt,mach_numnt,massraten1dt,massraten50dt,massraten100dt,massraten200dt,massraten500dt,massraten700dt,massraten1400dt] = fun_nonconser(L,n,x,dx,gamma,k);
execution_time_nonconser=toc
tic;
 [rhoc,tc,vc,Pc,mach_numc,rhoct,Pct,tct,mach_numct,massrate1dt,massrate50dt,massrate100dt,massrate200dt,massrate500dt,massrate700dt,massrate1400dt] = fun_conser(L,n,x,dx,gamma,k);
 execution_time_conser=toc

  figure(3)

  %comparison of actual values of conservative vs nonconservative
   subplot(4,1,1)
   plot(x,mach_numc)
   hold on
   plot(x,mach_num,\'*\')
   hold on
   xlabel(\'length of nozzle\')
   ylabel(\'mach number\')
   legend(\'conservative\',\'nonconservative\')
   title(\'comparison of actual values mach,P,rho,temp\')
   subplot(4,1,2)
   plot(x,Pc)
   hold on
   plot(x,P,\'*\')
   hold on
   xlabel(\'length of nozzle\')
   ylabel(\'pressure ratio\')
   legend(\'conservative\',\'nonconservative\')
   subplot(4,1,3)  
   plot(x,rhoc)
   hold on
   plot(x,rho,\'*\')
   hold on
   xlabel(\'length of nozzle\')
   ylabel(\'density ratio\')
    legend(\'conservative\',\'nonconservative\')
   subplot(4,1,4)
   plot(x,tc)
   hold on
   plot(x,t,\'*\')
   hold on
   xlabel(\'length of nozzle\')
   ylabel(\'temp ratio\')
    legend(\'conservative\',\'nonconservative\')

figure(4)
   %comparison of throat values along the nozzle
   subplot(4,1,1)
   plot(1:1400,mach_numct)
   hold on
    plot(1:1400,mach_numnt,\'r\')
   hold on
   xlabel(\'steps\')
   ylabel(\'mach number \')
   legend(\'conservative\',\'nonconservative\')
   title(\'comparison of throat values along the nozzle\')
   subplot(4,1,2)
   plot(1:1400,Pct)
   hold on
    plot(1:1400,Pnt,\'r\')
   hold on
   xlabel(\'steps\')
   ylabel(\'pressure ratio\')
   legend(\'conservative\',\'nonconservative\')
   subplot(4,1,3)  
    plot(1:1400,rhoct)
    hold on
    plot(1:1400,rhont,\'r\')
    hold on
   xlabel(\'steps\')
   ylabel(\'density ratio\')
    legend(\'conservative\',\'nonconservative\')
   subplot(4,1,4)
   plot(1:1400,tct)
   hold on
   plot(1:1400,tnt,\'r\')
   hold on
   xlabel(\'steps\')
   ylabel(\'temp ratio\')
    legend(\'conservative\',\'nonconservative\')
    
    
    
 figure(5)
    plot(x,massrate1dt)
    hold on
    xlabel(\'length of nozzle\')
    ylabel(\'mass flowrate\')
    title(\'mass flow rate at different time steps of conservative\')
    
    plot(x,massrate50dt)
    hold on
    xlabel(\'length of nozzle\')
    ylabel(\'mass flowrate\')
    
    plot(x,massrate100dt)
    hold on
    xlabel(\'length of nozzle\')
    ylabel(\'mass flowrate\')
    
    plot(x,massrate200dt)
    hold on
    xlabel(\'length of nozzle\')
    ylabel(\'mass flowrate\')
    
    plot(x,massrate500dt)
    hold on
    xlabel(\'length of nozzle\')
    ylabel(\'mass flowrate\')
    
    plot(x,massrate700dt)
    hold on
    xlabel(\'length of nozzle\')
    ylabel(\'mass flowrate\')
    
    plot(x,massrate1400dt)
    xlabel(\'length of nozzle\')
    ylabel(\'mass flowrate\')
    hold on
    hold off
    legend(\'1dt\',\'50dt\',\'100dt\',\'200dt\',\'500dt\',\'700dt\',\'1400dt\')

      
        
    figure(6)
    plot(x,massraten1dt)
    hold on
    xlabel(\'length of nozzle\')
    ylabel(\'mass flowrate\')
    title(\'mass flow rate at different time steps of non conservative\')
    
    plot(x,massraten50dt)
    hold on
    xlabel(\'length of nozzle\')
    ylabel(\'mass flowrate\')
    
    plot(x,massraten100dt)
    hold on
    xlabel(\'length of nozzle\')
    ylabel(\'mass flowrate\')
    
    plot(x,massraten200dt)
    hold on
    xlabel(\'length of nozzle\')
    ylabel(\'mass flowrate\')
    
    plot(x,massraten500dt)
    hold on
    xlabel(\'length of nozzle\')
    ylabel(\'mass flowrate\')
    
    plot(x,massraten700dt)
    hold on
    xlabel(\'length of nozzle\')
    ylabel(\'mass flowrate\')
    
    plot(x,massraten1400dt)
    xlabel(\'length of nozzle\')
    ylabel(\'mass flowrate\')
    hold on
    hold off
    legend(\'1dt\',\'50dt\',\'100dt\',\'200dt\',\'500dt\',\'700dt\',\'1400dt\')

Comparison between conservative forms and non-conservative forms :

Conclusion and result:

PLOTS:

1)Comparison of actual values as Mach number, Pressure ratio, temperature ratio, density ratio with length of nozzle:

2)Comparison of throat values along with the nozzle:

From the above plots, we can say that the timewise change of variables at the throat is much less for the conservative form than the nonconservative form. In both method at initial steps, changes are visible.

3)Mass flow rate at different time steps of conservative:

4)Mass flow rate at different time steps of nonconservative:

From the above plots, we can say that the variation in mass flow rates is less in conservative form. Hence the solution of conservative is accurate than nonconservative form.

  • After comparing both the conservative and nonconservative forms on 1d supersonic nozzle flow we get that most of the results are similar.
  • The conservation form yielded a slightly better mass flow distribution.

Comparing execution time(without any plots) of both the method we found:

The execution time of nonconservative: 0.1102

The execution time of conservative: 0.1705

Hence nonconservative is faster. But we also need to consider that we get accurate results in conservative method. After comparing both user can decide which method is required at a particular situation.

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.