SIMULATION OF 1D quasi super-sonic nozzle flow using maccormack methord

objective:

solve the 1D supersonic nozzle flow equations using the Macormack Method by implementing both the conservative and non-conservative forms of the governing equation.also:

• Figuring out the minimum number of cycles for which the simulation should be run in order for convergence
• Compararing the normalized mass flow rate between the conservative forms and non-conservative forms - making suitable conclusions
• Which method is faster

OUTPUTS AND OBSERVATIONS , PLOTS:

PLOT1: MASS FLOW RATE

conservative:

non-vonservative:

* we see that there are some oscillations at outflow and infloew boundaries of the non conservative form whereas in conservative form we see a constant value for mass flow with very litte variation.

analytical mass flow rate =rho*area*velocity

the conservative form is more close to the analytical solution of 0.579

PLOT2: MASS FLOW RATE per unit area

conservative

:

non-vonservative:

the plots for mass flow rate per unit area is maximum for at the throat for both conservative and non conservative cases ie 1.5m

PLOT3: VARIATION OF PROPERTIES ALONG THE LENGTH OF THE NOZZLE

conservative:

non-vonservative:

PLOT4: MACH NUMBER

conservative:

non-vonservative:

mach number value is 1 at the throat of the nozzle in both cases as we go from zero to 1.5 m the mach no increases upto 1 at x=1.5 mach no becomes 1 at the throat from x=1.5 to 3 m the mach no >1 and keeps on increasing

PLOT5: CONVERGENCE AT THE THROAT

conservative:

NOTE:the following lines represent

YELLOW :MACH NO

RED:TEMPERATURE(T/T*)

BLUE:DENSITY(RHO/RHO*)

GREEN:PRESSURE(P/P*)

non-vonservative:

NOTE:the following lines represent

YELLOW :MACH NO

RED:TEMPERATURE(T/T*)

BLUE:DENSITY(RHO/RHO*)

GREEN:PRESSURE(P/P*)

SIMULATION TIME TAKEN

the minimum time taken to converge by the two methords are

1> non conservative: 5.5

at iteration no 276 the error dec bolow tolerance level ie 10e-4;
hence time for reaching convergence is equal to 276*0.0201=5.5476
seconds approximately

2>conservative:3.7

at iteration no 182 the error dec bolow tolerance level ie 10e-4;
hence time for reaching convergence is equal to 182*0.0201=3.7
seconds approximately

hence we conclude that the time taken by conservative form takes less time as compared to the non conservative form of the governinng equation

and hence saves us a little bit of time whie processing

MINIMUM NO OF CYCLES FOR CONVERGENCE:

this can be determined by courants no

dt=C*dx/(a+v);

where a=sqrt(t);

v=velocity

dx=spacing between grid points

for the given simulation error tolerance is taken as 1e-4

gris spacing =0.1

no of grid points =31

calculating min no of cycles :

non conservative-->942

conservative-->819

GRID DEPENDENCE:

6 cases were concidered and mass flow rate was seen as a parameter for grid dependence from the given table was formed in ms word:

Grid points	conservative form	mass flow rate	non conservative form	mass flow rate
31	0.5853	819	0.5837	940
61	0.5839	1611	0.5825	1716
91	0.5824	2298	0.5809	2533
151	0.5821	3665	0.5815	3819
201	0.5821	4895	0.5819	7399
301	0.5822	6575	0.5819	8419
analytical mass flow rate	0.579

from the above table

we see that 31 grid points are more that enough for a feasable solution.

we see that it becomes grid size independent after 151 grid point and there is no major change by changing the grid size of both the forms.

code for conservative form:

 
clear all
close all
clc

n=31;
x=linspace(0,3,n);
dx=x(2)-x(1);
gamma=1.4;
iter1=linspace(0,1000,1);

tol=10e-4;

for i=1:n
if x(i)<=0.5
rho(i)=1;
t(i)=1;  %t is temperature
else if (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);
    else 
        rho(i)=0.634-0.3879*(x(i)-1.5);
        t(i)=0.833-0.3507*(x(i)-1.5);
    end
end
end

a=1+2.2*(x-1.5).^2;

%velocity
v=0.59./(rho.*a);

c=0.5 ; %couran number
del_t=c*dx./(sqrt(t)+v);
dt=min(del_t)%comes out to be 0.0102
nt=500;

%macormac methord
u1=rho.*a;
u2=rho.*a.*v;
u3=rho.*a.*((t/(gamma-1))+((gamma/2)*(v.^2)));

for k=1:nt
    rho_old=rho;
    v_old=v;
    t_old=t;
u1_old=u1;
u2_old=u2;
u3_old=u3;
  %  flux variables
  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 i=2:n-1
             dadx=((a(i+1))-a(i))/dx;
             j2=(1/gamma)*rho(i)*t(i)*dadx;
             
             %predictor step
             du1dt_pr(i)=-(f1(i+1)-f1(i))/dx;
             du2dt_pr(i)=-((f2(i+1)-f2(i))/dx)+j2;
             du3dt_pr(i)=-((f3(i+1)-f3(i))/dx);
             
             %updating value
             u1(i)=u1_old(i)+du1dt_pr(i)*dt;
             u2(i)=u2_old(i)+ du2dt_pr(i)*dt;
             u3(i)=du3dt_pr(i)*dt+u3_old(i);
        end
    %corrector step
    rho=u1./a;
    v=u2./u1;
    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 i=2:n-1
             
             dadx=(a(i)-a(i-1))/dx; 
             j2=(1/gamma)*rho(i)*t(i)*dadx;
             
             
             du1dt_co(i)=-(f1(i)-f1(i-1))/dx;
             du2dt_co(i)=-((f2(i)-f2(i-1))/dx)+j2;
             du3dt_co(i)=-((f3(i)-f3(i-1))/dx);
    end         
     
       
        %averaging the result
     
        du1dt_avg=0.5.*(du1dt_pr+du1dt_co);
        du2dt_avg=0.5.*(du2dt_pr+du2dt_co);
        du3dt_avg=0.5.*(du3dt_pr+du3dt_co);
         
        %final values of u1 u2 u3 for perticular time step
        for i=2:n-1
            u1(i)=du1dt_avg(i)*dt+u1_old(i);
            u2(i)=du2dt_avg(i)*dt+u2_old(i);
            u3(i)=du3dt_avg(i)*dt+u3_old(i); 
        end
        
        %boundry conditions
        %inlet
        u2(1)=2*u2(2)-u2(3);
        u3(1)=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);
        
        %final values of the velocity pressure temperate etc
        rho=u1./a;
        v=u2./u1;
        t=(gamma-1)*((u3./u1)-((gamma/2)*v.^2));
        
        pressure=rho.*t;
        machno=v./(t.^0.5);
        massflowr=v.*a.*rho;
        mfla=v.*rho;
        
    error1=max(abs(v-v_old))    
    error2=max(abs(rho_old-rho));
    error3=max(abs(t_old-t));
    error=max([error1;error2;error3]);
    
      if error

CODE FOR NOT CONSERVATIVE FORM:

clear 
close all
clc

%inputs
n=31;
x=linspace(0,3,n);
dx=x(2)-x(1);
gamma=1.4;
iter1=linspace(0,1000,1);


%calculating initial profiles
rho=1-0.3146*x;
t=1-0.2314*x;  %t is temperature
v=(0.1+1.09*x).*t.^0.5;
a=1+2.2*(x-1.5).^2;

%calculating dt using the couran no approach 
c=0.5 ; %couran number
del_t=c*dx./(sqrt(t)+v);
dt=min(del_t)

rho_ini=rho;
t_ini=t;
a_ini=a;
v_ini=v;
nt=1000;
error=9e9;
tol=10e-4;


for k=1:nt
    
    rho_old=rho;
    v_old=v;
    t_old=t;
    %predictor methord writing the continuity momentum and energy equation
    %from theory
    
    %continuity equation:
    %drhodt=-rho*dvdx-rho*v*(dloga)dx-v*()drhodx)
    
    %momentum equation 
    %dvdy=-v*dvdx-(1/gamma)*dtdx+t*drhodx/rho)
    
    %energy equation
    %dtdt=-v*dtdx-(gamma-1)*t(dvdx+v*dlogadx)
    
    for j=2:n-1
        dvdx=(v(j+1)-v(j))/dx;
        dlogadx=(log(a(j+1))-log(a(j)))/dx;
        drhodx=(rho(j+1)-rho(j))/dx;
        dtdx=(t(j+1)-t(j))/dx;
        drhodt_pr(j)=-rho(j)*dvdx-rho(j)*v(j)*dlogadx-v(j)*drhodx;
        
        dvdt_pr(j)=-v(j)*dvdx-(1/gamma)*(dtdx+t(j)*drhodx/rho(j));

        dtdt_pr(j)=-v(j)*dtdx-(gamma-1)*t(j)*(dvdx+v(j)*dlogadx); 
        
        %solution update 
        rho(j)=rho_old(j)+drhodt_pr(j)*dt;
        v(j)=v_old(j)+dvdt_pr(j)*dt;
        t(j)=t_old(j)+dtdt_pr(j)*dt;
        
    end
    
    %corrector methord
    for j=2:n-1
        dvdx=(v(j)-v(j-1))/dx;
        dlogadx=(log(a(j))-log(a(j-1)))/dx;
        drhodx=(rho(j)-rho(j-1))/dx;
        dtdx=(t(j)-t(j-1))/dx;
        
        drhodt_co(j)=-rho(j)*dvdx-rho(j)*v(j)*dlogadx-v(j)*drhodx;
        
        dvdt_co(j)=-v(j)*dvdx-(1/gamma)*(dtdx+t(j)*drhodx/rho(j));

        dtdt_co(j)=-v(j)*dtdx-(gamma-1)*t(j)*(dvdx+v(j)*dlogadx);
    end
    
    %calculating average derivative 
    drhodt_av=0.5*(drhodt_pr+drhodt_co);
    dtdt_av=0.5*(dtdt_pr+dtdt_co);
    dvdt_av=0.5*(dvdt_pr+dvdt_co);
    
    %next time step
   
        for i=2:n-1
    rho(i)=rho_old(i)+drhodt_av(i)*dt;
    v(i)=v_old(i)+dvdt_av(i)*dt;
    t(i)=t_old(i)+dtdt_av(i)*dt;
      %mass flow rate =density*velocity*area
      mfl(i)=rho(i)*a(i)*v(i);
      
      %mass flow rate per unit area computed =density*velocity
      mfla(i)=rho(i)*v(i);
      
      %mach number computed
      mach(i)=v(i)/(sqrt(t(i)));
      
      %pressure =density *temperature
      pressure(i)=rho(i)*t(i);
      
     
      
        end
    %boundry conditions 
    %inflow
    v(1)=2*v(2)-v(3);
    rho(1)=1;
    t(1)=1;
    pressure(1)=2*pressure(2)-pressure(3);
    mach(1)=2*mach(2)-mach(3);
    mfl(1)=2*mfl(2)-mfl(3);
    mfla(1)=2*mfla(2)-mfla(3);
    
    %outflow
    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);
    pressure(n)=2*pressure(n-1)-pressure(n-2);
    mach(n)=2*mach(n-1)-mach(n-2);    
    mfl(n)=2*mfl(n-1)-mfl(n-2);
    mfla(n)=2*mfla(n-1)-mfla(n-2);
     
       
    title(sprintf(\'time step = %d\',k));
    error1=max(abs(v_old-v));
    error2=max(abs(rho_old-rho));
    error3=max(abs(t_old-t));
    error=max([error1;error2;error3]);
    
      if error