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

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

Introduction:-

nozzles are most commonly used in every fields in order to study the liquids and air properties, because those two elements plays an important role in the enginnering fields. However we have many properties and the beheaviour of liquids and air but we are still in position to find the exact solutions for the equation which may help to find out the answers to unknown questions.

The project which based on nozzle in an important project to understand the working of the fluid particles when it is travelling through a convergent-divergent nozzle.

In this project we are using the macormack method to solve the equations because this method is widely used by the engineers to solve the complex second order derivatives, so obviously the complexicity of the equation is reduced.

 

Objective:-

The main objective of  this project is to learn the properties of a fluid when it is travelling through a convergent-divergent nozzle. 

As mentioned in the question there are two types of flows namely -

1. conservative flow

2. Non-conservative flow

 

CONSERVATIVE FLOW:-

In this conservative flow, when a fluid is flowing through a pipe a small control volume is taken in some area inside that pipe and this control volume is fixed, so that the equation can be solved to find out the properties of the flow at that point, this in known as the conservative flow.

NON-CONSERVATIVE FLOW:-

In this Non-conservative flow, a small control volume is taken inside the pipe flow but with a small change, the control volume here is moving with the flow instead of being fixed as we see in the conservative flow.

 

Programming concepts:-

We use the maccormac method for solving the equation in which the concept is used for programming to give exact solution. We use to kinds of programming methods in which these two methods are compared to find out that which method produce the answer exactly and quickly solution. These two methods are going to solve the same equation in order to find out that which method is suitable for the solution.

Initially the grid points are choosen as minimum value and then the gird points are increased to see the change in the solution. At first the grid points are taken as 31 because this is the maximum value for steady state to attain. At this point the graphs are plotted and the results is given below.

Secondly the grid points are made double (61) to check wheather there is a drastic change or not. These graphs  are also plotted and shown below.

In a state Both provide the same plots after the steady state is achieved. Nearly 1000 time steps are followed by time marching procedure to provide a steady state. Now lets see the time wise variation along throat side.

 

Variation of properties along the throat section:-

  • The concept we use in this method is one dimensional method. The whole problem is assumed as one dimensional material in which a point is considered to read the properties over that point.
  • but in reality there is always a two dimensional concept which flows through the nozzle.
  • By taking all these constraints in the program we are provided with the results to find out the exact solution.
  • As we see that Non-conservative flow provides more accurate results for primitive variables.
  • But in conservative flow we use many different substitution for solving the same.
  • As a result the non-conservative method provides less no. of error when compared to the conservative solutions.

Now in order to compare the values we are assuming different grid size, as we double the grid size the time steps should also be doubled in order to provide the exact solution for both the conservative and non conservative solutions.

 

PROGRAM FOR CONSERVATIVE AND NON-CONSERVATIVE FORM FLOW:-

1. Matlab function for Conservative nozzle flow:-

function [rho_c,t_c,m_c,mach_c,v_c,p_c,rho_c_t,t_c_t,m_c_t,mach_c_t,v_c_t,p_c_t]=conservative_form_flow(n,nt,trt)
%inputs
x=linspace(0,3,n);
dx=x(2)-x(1);
gamma=1.4;

%calculate initial profiles
for i= 1:n
    if (x(i)>=0 && x(i)<0.5)
        rho_c(i)=1;
        t_c(i)=1;
    elseif (x(i)>=0.5 && x(i)<1.5)
        rho_c(i)=1-(0.366*(x(i)-0.5));
        t_c(i)=1-0.167*(x(i)-0.5);
    elseif (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

%area
a=1 + 2.2*((x-1.5).^2);
c=0.5;
%velocity
v_c=0.59./(rho_c.*a);

%number of tiem steps


%solution vectors
p_c=rho_c.*t_c;
u_1=rho_c.*a;
u_2=rho_c.*a.*v_c;
u_3=u_1.*((t_c./(gamma-1))+((gamma/2).*(v_c.^2)));

%outer time loop
for k=1:nt
    dt=min (c.*(dx./(sqrt(t_c)+v_c)));
    %predictor method
    u_1old=u_1;
    u_2old=u_2;
    u_3old=u_3;
    %now in this conservative form a new term called flux term was
    %introduced and it should be also initialised for the prolem
    f_1=u_2;
    term=(((gamma-1)/gamma)*(u_3-((gamma/2)*((u_2.^2)./u_1))));
    f_2=((u_2.^2)./u_1) + term;
    f_3=((gamma*u_2.*u_3)./u_1)-(gamma*(gamma-1).*(u_2.^3))./(2*(u_1.^2));
    
    for j=2:n-1
        df_1dx(j)=(f_1(j+1)-f_1(j))/dx;
        df_2dx(j)=(f_2(j+1)-f_2(j))/dx;
        df_3dx(j)=(f_3(j+1)-f_3(j))/dx;
        %dlogadx(j)=(log(a(j+1))-log(a(j)))/dx;
        dadx(j)=(a(j+1)-a(j))/dx;
        j_2=(1/gamma)*rho_c(j).*t_c(j).*dadx(j);
        
        
        %continuity equation 
        du_1dt_p(j)= -df_1dx(j);
        
        %momentum equation 
        du_2dt_p(j)= -df_2dx(j) + j_2;
        
        %energy equation
        du_3dt_p(j)= -df_3dx(j);
        
        
        %solution update
        u_1(j)=u_1(j) + du_1dt_p(j)*dt;
        u_2(j)=u_2(j) + du_2dt_p(j)*dt;
        u_3(j)=u_3(j) + du_3dt_p(j)*dt;
        
            rho_c(j) = u_1(j)./a(j);
    v_c(j) = u_2(j)./u_1(j);
    t_c(j) = (gamma-1)*((u_3(j)./u_1(j))-((gamma*(v_c(j).^2))/2));
        
    end
    
%     rho_c = u_1./a;
%     v_c = u_2./u_1;
%     t_c = (gamma-1)*((u_3./u_1)-((gamma*v_c.^2)/2));
    
    f_1=u_2;
    term_1=(((gamma-1)/gamma)*(u_3 - (gamma*(u_2.^2))./(2*u_1)));
    f_2=((u_2.^2)./u_1)+term_1;
    f_3=((gamma*u_2.*u_3)./u_1)-(gamma*(gamma-1)*(u_2.^3))./(2*(u_1.^2));
    
    %corrector step
       
    for j=2:n-1
        df_1dx(j)=(f_1(j)-f_1(j-1))/dx;
        df_2dx(j)=(f_2(j)-f_2(j-1))/dx;
        df_3dx(j)=(f_3(j)-f_3(j-1))/dx;
%         dlogadx(j)=(log(a(j))-log(a(j-1)))/dx;
        dadx(j)=(a(j)-a(j-1))/dx;
        j_3(j)=(1/gamma)*rho_c(j)*t_c(j)*dadx(j);
        
        %continuity equation 
        du_1dt_c(j)=-df_1dx(j);
        
        %momentum equation 
        du_2dt_c(j)=-df_2dx(j)+j_3(j);
        
        %energy equation
        du_3dt_c(j)=-df_3dx(j);
    end
   %compute the avarage time derivative
   du_1dt=0.5*(du_1dt_p+du_1dt_c);
   du_2dt=0.5*(du_2dt_p+du_2dt_c);
   du_3dt=0.5*(du_3dt_p+du_3dt_c);
  
   for i=2:n-1
       %final solution update
        u_1(i)=u_1old(i)+du_1dt(i)*dt;
        u_2(i)=u_2old(i)+du_2dt(i)*dt;
        u_3(i)=u_3old(i)+du_3dt(i)*dt;
   end
   %applying the boundary conditions
   %inlet BC
   u_1(1)=rho_c(1)*a(1);
   u_2(1)=2*u_2(2)-u_2(3);
   v_c(1)=u_2(1)/u_1(1);
   u_3(1)=u_1(1)*((t_c(1)/(gamma-1))+((gamma/2)*(v_c(1)^2)));
   
   %outlet BC
   u_1(n)=2*u_1(n-1)-u_1(n-2);
   u_2(n)=2*u_2(n-1)-u_2(n-2);
   u_3(n)=2*u_3(n-1)-u_3(n-2);
   
   rho_c = u_1./a;
   v_c = u_2./u_1;
   t_c = (gamma-1)*((u_3./u_1)-((gamma*v_c.^2)/2));
 
   %calculating 
   mach_c=v_c ./ (t_c.^0.5);
   m_c=rho_c .* v_c .* a;
   p_c=rho_c .*t_c;
   
   %at throat
   rho_c_t(k)=rho_c(trt);
   p_c_t(k)=p_c(trt);
   t_c_t(k)=t_c(trt);
   mach_c_t(k)=mach_c(trt);
   v_c_t(k)=v_c(trt);
   m_c_t(k)=m_c(trt);
   
   %now plotting the mass flow rate at different time steps
   figure(1)
   if k==50
       plot(x,m_c,'r','linewidth',1.5)
       hold on
   elseif k==100
        plot(x,m_c,'b','linewidth',1.5)
        hold on
   elseif k==200
        plot(x,m_c,'g','linewidth',1.5)
        hold on
   elseif k==400
        plot(x,m_c,'y','linewidth',1.5)
        hold on
   elseif k==500
        plot(x,m_c,'c','linewidth',1.5)
        hold on
   elseif k==700
        plot(x,m_c,'k','linewidth',1.5)
        hold on
        
        grid on
        xlabel('Non-dimensional distance');
        ylabel('Non-dimensional mass flow');
        title('Conservative form of nozzle Flow')
        legend({'50 time steps','100 time steps','200 time steps','400 time steps','500 time steps','700 time steps'});
   end
   
end
end

 

 

2. Matlab function for Non-conservative nozzle flow:-

function [rho_n,t_n,m_n,mach_n,v_n,p_n,rho_n_t,t_n_t,m_n_t,mach_n_t,v_n_t,p_n_t]=non_conservative_form(n,nt,trt)
%inputs
x=linspace(0,3,n);
dx=x(2)-x(1);
gamma=1.4;

%calculate initial profiles
rho_n=1-0.3146*x;
t_n=1-0.2314*x;  %t=temperature
v_n=(0.1+1.09*x).*t_n.^0.5;

%area
a=1+2.2*(x-1.5).^2;
c=0.5;
%number of tiem steps
dt=min(c.*(dx./((t_n.^0.5)+v_n)));

%outer time loop
for k=1:nt
    
    %predictor method
    rho_old=rho_n;
    v_old=v_n;
    t_old=t_n;
  
    for j=2:n-1
        %drho_dt_p = -(rho(j))*((v(j+1)-v(j))/dx)-rho(j)*v(j)*((log(a(j+1))-log(a(j)))/dx)-v(j)*((rho(j+1)-rho(j))/dx);
        dvdx=(v_n(j+1)-v_n(j))/dx;
        drhodx=(rho_n(j+1)-rho_n(j))/dx;
        dtdx=(t_n(j+1)-t_n(j))/dx;
        dlog_a_dx=(log(a(j+1))-log(a(j)))/dx;
        %continuity equation
        drho_dt_p(j)=-rho_n(j)*dvdx -rho_n(j)*v_n(j)*dlog_a_dx -v_n(j)*drhodx;
        
        %momentum equation
        dvdt_p(j) = -v_n(j)*dvdx -(1/gamma)*(dtdx + (t_n(j)/rho_n(j))*drhodx);
        
        %energy equation
        dtdt_p(j)=-v_n(j)*dtdx -(gamma-1)*t_n(j)*(dvdx + v_n(j)*dlog_a_dx);

        %solution update
         rho_n(j)=rho_n(j)+drho_dt_p(j)*dt;
         v_n(j)=v_n(j)+dvdt_p(j)*dt;
         t_n(j)=t_n(j)+dtdt_p(j)*dt;
    end
    
    %corrector method

    for j=2:n-1
        %drho_dt_p = -(rho(j))*((v(j)-v(j-1))/dx)-rho(j)*v(j)*((log(a(j))-log(a(j-1)))/dx)-v(j)*((rho(j)-rho(j-1))/dx);
        dvdx=(v_n(j)-v_n(j-1))/dx;
        drhodx=(rho_n(j)-rho_n(j-1))/dx;
        dtdx=(t_n(j)-t_n(j-1))/dx;
        dlog_a_dx=(log(a(j))-log(a(j-1)))/dx;
        %continuity equation
        drho_dt_c(j)=-rho_n(j)*dvdx-rho_n(j)*v_n(j)*dlog_a_dx-v_n(j)*drhodx;
        
        %momentum equation
        dvdt_c(j) = -v_n(j)*dvdx-(1/gamma)*(dtdx+(t_n(j)/rho_n(j))*drhodx);
        
        %energy equation
        dtdt_c(j)=-v_n(j)*dtdx-(gamma-1)*t_n(j)*(dvdx+v_n(j)*dlog_a_dx);

    end
   %compute the avarage time derivative
   drhodt=0.5*(drho_dt_p+drho_dt_c);
   dvdt=0.5*(dvdt_p+dvdt_c);
   dtdt=0.5*(dtdt_p+dtdt_c);
  
   for i=2:n-1
       %final solution update
        rho_n(i)=rho_old(i)+drhodt(i)*dt;
        v_n(i)=v_old(i)+dvdt(i)*dt;
        t_n(i)=t_old(i)+dtdt(i)*dt;
   end
   %applying the boundary conditions
   %inlet BC
   v_n(1)=2*v_n(2)-v_n(3);
   %outlet BC
   rho_n(n)=2*rho_n(n-1)-rho_n(n-2);
   v_n(n)=2*v_n(n-1)-v_n(n-2);
   t_n(n)=2*t_n(n-1)-t_n(n-2);
   
      %calculating 
   mach_n=v_n./(t_n.^0.5);
   m_n=rho_n.*v_n.*a;
   p_n=rho_n.*t_n;
   
   %at throat
   rho_n_t(k)=rho_n(trt);
   p_n_t(k)=p_n(trt);
   t_n_t(k)=t_n(trt);
   mach_n_t(k)=mach_n(trt);
   v_n_t(k)=v_n(trt);
   m_n_t(k)=m_n(trt);
   
   %now plotting the mass flow rate at different time steps
   figure(2)
   if k==50
       plot(x,m_n,'r','linewidth',1.5)
       hold on
   elseif k==100
        plot(x,m_n,'b','linewidth',1.5)
        hold on
   elseif k==200
        plot(x,m_n,'g','linewidth',1.5)
        hold on
   elseif k==400
        plot(x,m_n,'y','linewidth',1.5)
        hold on
   elseif k==500
        plot(x,m_n,'c','linewidth',1.5)
        hold on
   elseif k==700
        plot(x,m_n,'k','linewidth',1.5)
             
        grid on
        xlabel('Non-dimensional distance');
        ylabel('Non-dimensional mass flow');
        title('Non-Conservative form of nozzle Flow')
        legend({'50 time steps','100 time steps','200 time steps','400 time steps','500 time steps','700 time steps'});
   end
end  
end

 

3. Main program for nozzle flow:-

clear all
close all
clc

%now providing initial condition
n=31;
x=linspace(0,3,n);

%now providing the node values
trt=((n-1)/2)+1;

%number of timesteps
nt=700;

tic
%now providing the function of conservative form
[rho_c,t_c,m_c,mach_c,v_c,p_c,rho_c_t,t_c_t,m_c_t,mach_c_t,v_c_t,p_c_t]=conservative_form_flow(n,nt,trt);
time_c=toc;
tic
%now providing the function for the non conservative programe 
[rho_n,t_n,m_n,mach_n,v_n,p_n,rho_n_t,t_n_t,m_n_t,mach_n_t,v_n_t,p_n_t]=non_conservative_form(n,nt,trt);
time_n=toc;

%plotting 
 figure(3)
 subplot(4,1,1)
 plot(linspace(1,nt,nt),rho_c_t);
 hold on 
 plot(linspace(1,nt,nt),rho_n_t,'r');
 %hold on 
 ylabel('rho')
 legend({'Conservative form','Non conservative form'})
 axis([0 1500 0 1.5]);
 title('time variation at the throat')
 
 
 subplot(4,1,2)
 plot(linspace(1,nt,nt),t_c_t);
 hold on 
 plot(linspace(1,nt,nt),t_n_t,'r');
 %hold on 
 ylabel('t-temperature')
 legend({'Conservative form','Non conservative form'})
 axis([0 1500 0 1.5]);
 
 
 subplot(4,1,3)
 plot(linspace(1,nt,nt),p_c_t);
 hold on 
 plot(linspace(1,nt,nt),p_n_t,'r');
 %hold on 
 ylabel('p-pressure')
 legend({'Conservative form','Non conservative form'})
 axis([0 1500 0 1.5]);
 
 
 subplot(4,1,4)
 plot(linspace(1,nt,nt),mach_c_t);
 hold on 
 plot(linspace(1,nt,nt),mach_n_t,'r');
 %hold on 
 ylabel('mach number')
 legend({'Conservative form','Non conservative form'})
 axis([0 1500 0 1.5]);
 
 
 
 figure(4)
 subplot(4,1,1)
 plot(x,p_c);
 hold on
 plot(x,p_n,'--','color','k')
 grid on
 ylabel('pressure')
 legend({'conservative form','non-conservative form'})
 
 subplot(4,1,2)
 plot(x,rho_c);
 hold on
 plot(x,rho_n,'--','color','k')
 grid on
 ylabel('density')
 legend({'conservative form','non-conservative form'})
 
 subplot(4,1,3)
 plot(x,t_c);
 hold on
 plot(x,t_n,'--','color','k')
 grid on
 ylabel('temperature')
 legend({'conservative form','non-conservative form'})
 
 subplot(4,1,4)
 plot(x,mach_c);
 hold on
 plot(x,mach_n,'--','color','k')
 grid on
 ylabel('mach no.')
 legend({'conservative form','non-conservative form'})
 
 
 %mass plow rate comperission 
 figure(5)
 plot(x,m_c,'b')
 hold on
 plot(x,m_n,'y')
 hold on
 ylabel('non dimensional mass flow rate')
 xlabel('non dimensional distance along x-axis')
 legend({'conservative','non-conservative'})
 grid on

 

Now lets see the  output graphs of different parameters :-

1. Mass flow rate for conservative nozzle flow:-

 

2. Mass flow rate for Non-conservative nozzle flow:-

 

3. variation of different properties with time variation:- 

 

4. variation of different properties  alog x axis :-

 

5. Mass flow rate comparision:-

 

Now as we see the above graphs clearly depicts the conservative and Non-conservative flow graphs results. These graphs are plotted for 31 grid points.

As we can see the conservative mass flow chart provide the more stable results when copared to the Non-conservative mass flow, because the mass flow is one of the dependent variables in the conservative flow.

So as a result, the conservative flow provides more accurate results for flux variables. So they are mostly used for conservative flow.

 

Now as we see the grid points in the above mentioned nozzle flow programs is 31, it is because about this grid size steady state has been occured for both conservative and non-conservative flow equations.

Now as we know the maximum grid points, lets check out the grid-independence for the conservative and non-conservative flow programs.

Grid independence is known as changing the grid points to check wheather the conservative and non-conservative flow may provide the same results same as before.

Now the grid points are doubled as 61, and time steps values are also increased to obtain the elapsed non-dimentional time for both the flow maybe the same. By increasing the time steps we can carry the grid independence by comparing the graphs obtained for both the cases.

 

Now let's see the Output graphs for increased gird size:-

1. mass flow rate for conservative nozzle flow:-

 

2.mass flow rate for Non-conservative nozzle flow:-

 

3. variation of different properties with time:-

 

4. variation of different properties along with x axis:-

 

5. mass flow rate comparision:-

 

Now as we can see the difference in thee graphs, we can notice that for n=61 the steady state is obtained after 1000 time steps.

 

Now which method is faster:-

We have made two different solutions for the equations, now we need to find that which method is faster than one another

For n=31:-

conservation time=1.8171

Non-conservation time=1.0713

 

For n=61:-

conservation time=2.4409

Non-conservation time=1.4780

 

So from these above data we came to a conclusion that conservation form requires more computational time than Non-conservative form.


Projects by Saurabh sharma

Challenge on Combustion
Saurabh sharma · 2020-02-16 18:45:07

  Challenge on Combustion Objective:- Perform a combustion simulation on the combuster model and plot the variation of the mass fraction of the different species in the simulation using line probes at different locations of the combuster. You need to plot Read more

 

Cyclone Separator Challenge

 

Objective:-

We have to perform an analysis on the cyclone seperator model for  different number of inectioned particles.

1. The four imperical models use Read more

  Gate Valve Parametric Study A gate valve, also known as a sluice valve, is a valve that opens by lifting a barrier(gate) out of the path of the fluid.Gate valves are used to shut off the flow of liquids rather than for flow regulation. Gate valve require Read more

Conjugate Heat Transfer Analysis on a graphics card. Description of a conjugate heat transfer analysis:- Conjugate heat transfer represents the combination of heat transfer in liquids and solids simulteneously. The mode of heat transfer in solids is conduction and con Read more

    Rayleigh Taylor Instability Challenge SOLUTION-1:- Models used for Rayleigh-Taylor instability:- Early-time growth Potential flow Buoyancy-drag balance Energy budget   1. Early-time growth:- It was noted in the taylor (1950) that t Read more

Exhaust Port Challenge
Saurabh sharma · 2019-12-19 08:14:22

              Exhaust Port Challenge SOLUTION1:- Why and Where a CHT analysis is used? Conjugate Heat Transfer analysis:- Conjugate heat transfer analysis is generally used when there is a teperature variation during heating and co Read more

Ahmed Body Challenge
Saurabh sharma · 2019-12-17 16:08:31

              Ahmed Body Challenge  SOLUTION1: Ahmed body:- The Ahmed body was originally described By S.R. Ahmed in 1984. Ahmed body is a simplified car body. The shape provides a  model to study geometric effects Read more

    Simulation of Flow through a pipe in                                OpenFoam Objective of the project:- 1. Calculation of flow variables. 2. Generation of the blockMeshDict Read more

             FVM:: Interpolation and Gradient Schemes literature review FVM-Flux limiters and interpolation schemes What is FVM? FVM is a finite volume method for representing and evaluating PDE\'s in the form of algebraic equations Read more


Loading...

The End