IN THIS PROJECT, I HAVE WORKED ON THE ONE DIMENSIONAL SUPERSONIC NOZZLE FLOW. IN THIS FLOW, I HAVE ASSUMED THE FLOW TO BE STEADY,ISENTROPIC. IN THIS NOZZLE, I HAVE CONSIDERED THE FLOW AT INLET SECTION COMES FROM RESERVOIR WHERE THE PRESSURE, TEMPERATURE ARE DENOTED AS `P_0` AND `T_0` RESPECTIVELY. THE CROSS-SECTIONAL AREA OF RESERVOIR IS INFINITE DUE TO WHICH THE VELOCITY BECOMES ZERO THERE.THE NOZZLE IS BASICALLY DIVIDED INTO TWO SECTIONS i.e CONVERGENT SECTION AND DIVERGENT SECTION. AS A RESULT, THE FLOW EXPANDS ISENTROPICALLY AND ATTAINS SUPERSONIC(M>1) SPEED AT EXIT OR DIVERGENT SECTION OF THE NOZZLE WHERE IT BECOMES SUBSONIC(M<1) AT THE CONVERGENT SECTION AND AT A POINT IT\'S SPEED  BECOMES EQUAL TO THE SONIC SPEED(M=1) WHICH IS KNOWN AS THROTTLE POINT(CENTRE POINT OF CONVERGENT AND DIVERGENT SECTION).

IT IS ALSO KNOWN AS THE QUASI-ONE DIMENSIONAL NOZZLE FLOW BECAUSE IN THE  FLOW, THE FLOW PROPERTIES ARE CONSIDERED AS THE FUNCTION OF DIRECTION OF THE FLOW.

THE GOVENING EQUATIONS USED IN THIS FLOW ARE:-

CONTINUITY EQUATION:-

                                  `rho_1*A_1+ rho_1*A_1*(v_1)^2+int_(A_1)^(A_2)rho*dA=rho_2*A_2+rho_2*A_2*(v_2)^2`

MOMENTUM EQUATION :-

                                   `h_1+(v_1)^2/2=h_2+(v_2)^2/2`

THE TECHNIQUE USED IN THIS NOZZLE FLOW FOR NUMERICAL SIMULATION SO AS TO FIND NATURE OF FLOW IS :-

MacCormack method:-

IN COMPUTATIONAL FLUID DYNAMICS, IT IS WIDELY USED IN DISCRETIZATION SCHEME FOR NUMERICAL SIMULATION OF HYPERBOLIC PARTIAL DIFFERENTIAL EQUATIONS.IT IS A SECOND ORDER DIFFERENCE METHOD.

The MacCormack method is a variation of the two-step Lax–Wendroff scheme but is much simpler in application. To illustrate the algorithm, consider the following first order hyperbolic equation

The application of MacCormack method to the above equation proceeds in two steps; a predictor step which is followed by a corrector step.

Predictor step: In the predictor step, a \"provisional\" value of  is estimated as follows

The above equation is obtained by replacing the spatial and temporal derivatives in the previous first order hyperbolic equation using forward differences.

Corrector step: In the corrector step, the predicted value  is corrected according to the equation

Note that the corrector step uses backward finite difference approximations for spatial derivative. Note also that the time-step used in the corrector step is in contrast to the  used in the predictor step.

Replacing the  term by the temporal average

to obtain the corrector step as

                                                                                                                        .

IN THIS PROJECT, I HAVE WORKED ON THE TWO FORMS OF NOZZLE FLOW EQUATION:-

1.) CONSERVATIVE FORM.

2.) NON-CONSERVATIVE FORM.

1.) NON-CONSERVATIVE FORM:- IN THIS FORM, WE HAVE CONSIDERED THE FLOW TO BE CONSEREVED i.e THE FLOW DOMAIN IS CONSTANT.

THE GOVERNING EQUATION IN THIS FORM IS:-

 CONTINUITY :- 

 `(del(rhoA))/(delt)+rhoA(delV)/(delx)+rhoV(delA)/(delx)+VA(delrho)/(delx)=0.`

 MOMENTUM :-

 `rho(delv)/(delt)+rhoV(delV)/(delx)=-R(rho(delT)/(delx)+T(delp)/(delx)).`

 ENERGY:-

 `rhoc_v(delT)/(delt)+rhoVc_v(delT)/(delx)=-rhoRT[(delV)/(delx)+V(del(lnA))/(delx)].`

HENCE, THE GOVERNING EQUATION IN THE FORM OF DIFFERENCE EQUATTION IS:-

`((delrho)/(delt))_i=-rho_i[(V_(i+1)-V_i)/(Deltax)]-rho_iV_i[lnA_(i+1)-lnA_i]/(Deltax)-V_i[rho_(i+1)-rho_i]/(Deltax).`

`((delV)/(delT))_i=-V_i[V_(i+1)-V_i]/(Deltax)-1/gamma[(T_(i+1)-T_i)/(Deltax)+T_i/rho_i((rho_(i+1)-rho_i)/(Deltax))].`

 `((delT)/(delt))_i=-V_i[(T_(i+1)-T_i)/(Deltax)]-(gamma-1)T_i[(V_(i+1)-V_i)/(Deltax)+V_i(lnA_(i+1)-lnA_i)/(Deltax)].`

2.) CONSERVATION FORM:- IN THIS FORM, WE HAVE CONSIDERED THE FLOW TO BE UNCONSEREVED i.e THE FLOW DOMAIN IS MOVING.

THE GOVERNING EQUATION IN THIS FORM IS:-

 `(del(rhoA))/(delt)+(del(rhoAV))/(delx)=0.`

`(del(rhoAV))/(delt)+(del[rhoAV+(1/gamma)rhoA])/(Deltax)=1/gammaP(delA)/dx.`

 `(del[rho(e^(gamma)/(gamma-1)+gamma/2V^2)A])/(Deltat)+(del[rho(e^gamma/(gamma-1)+gamma/2V^2)VA+rhoAV])/(Deltax)=0.`

NOW, IN ORDER TO RUN THE GOVERNING EQUATION, WE HAVE TO PROVIDE BOUNDARY CONDITION AND INTIAL CONDITION SO AS TO MAKE THE SIMULATION PROCESS MORE EASIER.

INTIAL CONDITION:-

 1.) NON-CONSERVATIVE:-

      FOR DENSITY,     `rho=1-0.3146x`

      FOR TEMPERATURE,  `T=1-0.2314X`

      FOR VELOCITY,   `V=(0.1+1.09x)T^(1/2)`

 2.) CONSERVATIVE :-

      `{(rho=1.0),(T=1.0):}}for0<=x<=0.5`

      `{(rho=1.0-0.366(x-0.5)),(T=1.0-0.167(x-0.5)):}}for0.5<=x<=1.5`

      `{(rho=0.634-0.3879(x-1.5)),(T=0.833-0.3507(x-1.5)):}}for1.5<=x<=3.5`

 NOZZLE SHAPE AND ITS CONDITION FOR BOTH :-

 `A=1+2.2(x-1.5)^2 , 0lexle3`

BOUNDARY CONDITION:-

SUBSONIC INFLOW BOUNDARY CONDITION:-

`V_1=2V_2-V_3`

`rho=1`

`T=1`

SUBSONIC OUTFLOW BOUNDARY CONDITION:-

`V_n=2*V_(n-1)-V_(n-2)`

`rho_n=2rho_(n-1)-rho_(n-2)`

`T_n=2T_(n-1)-T_(n-2)`

CALCULATION OF TIME STEP:- 

HERE TIME STEP ALSO PLAYS AN IMPORTANT ROLE IN BOTH CASE AS IT IMPACT IS HIGH ON THE CONVERGENCE OF VALUE . SO IN ORDER TO MAKE SIMULATION MORE CONVINIENT WE HAVE TO PROVIDE SUCH TIME STEP VALUE HERE WHICH WILL MAKE SIMULATION RUN EASIER AND WOULD MAKE US EASIER TO ATTAIN VALUE.

CRITERIA OF TIME STEP:-

`Deltat=C*((Deltax)/(a +v))`

where, \'a\' is known as as the speed of sound and speed of fluid is known as \'v\'.

\'C\' IS THE \'COURANT-NUMBER\'. WHERE IT SHOULD BE LESS THAN 1 FOR CONVERGENCE OF VALUE.

AFTER CALCULATING THE VALUE OF TIME STEP AT EACH GRID POINT OF NOZZLE, WE FIND OUT MINIMUM VALUE OF TIME STEP FROM IT WHICH WILL DO CONVERGENCE EASILY AT ALL GRID POINTS.

i.e `min(Deltat_1+Deltat_2+Deltat_3+....Deltat_n)`

MAIN CODE:-

clear all
close all 
clc

  
% Time steps:-
  nt=[100 400 700 1400];
  n=1;
  
 while(n|=0)
 
   n=input(\"Pick a number, any number for grid points for eg.31,61!\\n say 0 for no output\"); 
   x=linspace(0,3,n);
   dx=x(2)-x(1);
 
  for j=1:4
  [sim_tim_conv]=nozzle_flow_non_conservative(n,x,dx,nt(j));
  end
 
  for j=1:4
  [sim_tim_non_conv]=nozzle_flow_conservative(n,x,dx,nt(j));
  end 

  end

##### PLEASE INPUT VALUE IN COMMAND WINDOW  ####

FUNCTION DEFINED FOR CONSERVATIVE FORM:-

function [sim_tim]=nozzle_flow_conservative(n,x,dx,nt)

% calculation of Intial profile:-
  a=1+(2.2*(x-1.5).^2);
  
  for i=1:n
    if (x(i)>=0&&x(i)<=0.5)
       rho(i)=1;
       T(i)=1;
     end
    
     if (x(i)>=0.5&&x(i)<=1.5)
      rho(i)=1-(0.366*(x(i)-0.5));
      T(i)=1.0-(0.167*(x(i)-0.5));
     end
    
     if (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
      
  % Intial values:-
    v=0.59./(rho.*a);
    gamaa=1.4;
  
  % continuity equation:-
  
    U_1=rho.*a;
  
  % momentum equation:-
    U_2=rho.*a.*v;   
  
  % energy equation:- 
    U_3=rho.*((T./(gamaa-1))+((gamaa/2)*(v.^2))).*a;     

  % Time steps:-
    
    a_conv=sqrt(T);
    c=0.5;
    dt=(c*dx)./(a_conv+v);
    dt=min(dt);
    iter=0;
  
 tic()  
 % Time loop:-   
   for t=1:nt
  
    U_1_old=U_1;
    U_2_old=U_2;
    U_3_old=U_3;
    iter=iter+1;
  
 
 % predictor method:-
    
    % flux terms:-
    F_1=U_2;
   
    F_2=((U_2.^2)./U_1)+((gamaa-1)/gamaa)*(U_3-((gamaa/2)*((U_2.^2)./U_1)));
   
    F_3=gamaa*((U_2.*U_3)./U_1)-((gamaa*(gamaa-1))/2)*(U_2.^3./U_1.^2);
   


    for i=2:n-1
    
    % source term:-
    J_2(i)=(1/gamaa)*rho(i)*T(i)*(a(i+1)-a(i))/dx;
    
    
  % finding predicted values of solution terms using forward difference:-
    U_1_dt_p(i)= -(F_1(i+1)-F_1(i))/dx;
    
    U_2_dt_p(i)= -(F_2(i+1)-F_2(i))/dx+J_2(i);
    
    U_3_dt_p(i)= -(F_3(i+1)-F_3(i))/dx;
    
  % solution update:-
    U_1(i)=U_1(i)+ U_1_dt_p(i)*dt;
    
    U_2(i)=U_2(i)+ U_2_dt_p(i)*dt;
    
    U_3(i)=U_3(i)+ U_3_dt_p(i)*dt;
    
  end
  
   
  % updated predicted values of density and temperature:-
    rho=U_1./a;
    T=(gamaa-1)*(U_3./U_1-(gamaa/2)*((U_2.^2/U_1)));
    v=U_2./U_1; 
 
    % updated values of flux terms:-
    
    F_1=U_2;
   
    F_2=((U_2.^2)./U_1)+((gamaa-1)/gamaa)*(U_3-(gamaa/2)*(U_2.^2./U_1));
   
    F_3=gamaa*((U_2.*U_3)./U_1)-((gamaa*(gamaa-1))/2)*(U_2.^3./U_1.^2);
    
  
  % corrector step:-
    
    for i=2:n-1
     
    % updated value of source term:-
    
    J_2(i)=(1/gamaa)*rho(i)*T(i)*(a(i)-a(i-1))/dx;

  % finding solution terms using reaward differennce:-
    U_1_dt_c(i)= -(F_1(i)-F_1(i-1))/dx;
    
    U_2_dt_c(i)= -(F_2(i)-F_2(i-1))/dx +J_2(i);
    
    U_3_dt_c(i)= -(F_3(i)-F_3(i-1))/dx;   
   
    end

  
 % average time derivative is as folow:-
    
    U_1_dt_av= 0.5*(U_1_dt_p+U_1_dt_c);
    
    U_2_dt_av= 0.5*(U_2_dt_p+U_2_dt_c);
    
    U_3_dt_av= 0.5*(U_3_dt_p+U_3_dt_c);   
 
   % final solution update:-
   for j=2:n-1
      
    U_1(j)= U_1_old(j)+U_1_dt_av(j)*dt;
    
    U_2(j)= U_2_old(j)+U_2_dt_av(j)*dt;
    
    U_3(j)= U_3_old(j)+U_3_dt_av(j)*dt;
      
    end
  
## applying boundary conditions:-
  % inflow boundary condition:-
    U_2(1)=2*U_2(2)-U_2(3);
  
  % outflow boundary condition:-
   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);  
   
 end
 
 sim_tim=toc();


  % final value of density :-
    rho=U_1./a;
    
  % final value of temperature:-  
    T=(gamaa-1)*(U_3./U_1-(gamaa/2)*((U_2./U_1).^2)); 
    
  % final value of velocity:-
    v=U_2./U_1; 
    
  % final value of mass flow:-
    mass_flow_rate=rho.*a.*v;
    
  % final value of MACH number:-
    M=v./a_conv;
    
  % final value of pressure:-
    P=rho.*T;
    
  text=sprintf(\'CONSERVATIVE FLOW \\n NUMBER OF GRID POINTS USED IN NOZZLE LENGTH=%d,\\n NUMBER OF TIME STEPS=%d\',n,t);
    
  if(t==100) 
   figure(6)
   plot(x,mass_flow_rate,\'b\')
   axis([0 3 0.4 0.7]) 
   xlabel(\'length of nozzle (x/L)\')
   ylabel(\'mass flow rate(\\rho*a*v)/(\\rho_0*a_0*A)\') 
   title(text)
 end
 
   if(t==400) 
   figure(7)
   plot(x,mass_flow_rate,\'b\')
   axis([0 3 0.4 0.7]) 
   xlabel(\'length of nozzle (x/L)\')
   ylabel(\'mass flow rate(\\rho*a*v)/(\\rho_0*a_0*A)\') 
    title(text)
 end
 
   if(t==700) 
   figure(8)
   plot(x,mass_flow_rate,\'b\')
   axis([0 3 0.4 0.7]) 
   xlabel(\'length of nozzle (x/L)\')
   ylabel(\'mass flow rate(\\rho*a*v)/(\\rho_0*a_0*A)\') 
    title(text)
 end
 
   if(t==1400) 
   figure(9)
   plot(x,mass_flow_rate,\'b\')
   axis([0 3 0.4 0.7]) 
   xlabel(\'length of nozzle (x/L)\')
   ylabel(\'mass flow rate(\\rho*a*v)/(\\rho_0*a_0*A)\')
   title(text) 
   end

FUNCTION DEFINED FOR NON-CONSERVATIVE FORM:-

function [sim_tim]=nozzle_flow_non_conservative(n,x,dx,nt)

% calculation of Intial profile:-
  rho=1-0.3146*x;
  T=1-0.2314*x;
  v=(0.1+1.09*x).*T.^(1/2);
  a=1+2.2*(x-1.5).^2;
    
  gamaa=1.4;
 
% calculation of minimum time steps:- 
  a_conv=sqrt(T);
  c=1/2;
  dt=c*dx./(a_conv+v);
  dt=min(dt);
  
% intial iteration counter:-
  iter=0;
    
  tic()
  % Time loop:-
   for t=1:nt
     rho_old=rho;
     v_old=v;
     T_old=T;
     iter=iter+1;
          
  % predictor method:-
     
     for i=2:n-1
            
       dv_dx=(v(i+1)-v(i))/dx;
       drho_dx=(rho(i+1)-rho(i))/dx;
       da_dx=(log(a(i+1))-log(a(i)))/dx;
       dT_dx=(T(i+1)-T(i))/dx;
       
      % continuity equation
       drho_dt_p(i)=-rho(i)*dv_dx-rho(i)*v(i)*da_dx-v(i)*drho_dx;
       
      % momentum equation
       dv_dt_p(i)=-v(i)*dv_dx-(1/gamaa)*(dT_dx+(T(i)/rho(i))*drho_dx);
      
      % energy equation 
       dT_dt_p(i)=-v(i)*dT_dx-(gamaa-1)*T(i)*(dv_dx+v(i)*da_dx); 
       
       % Solution update
        v(i)=v(i)+dv_dt_p(i)*dt;
        rho(i)=rho(i)+drho_dt_p(i)*dt;
        T(i)=T(i)+dT_dt_p(i)*dt;
        
       end
      
   % corrector method:-  
   
       for i=2:n-1
       dv_dx=(v(i)-v(i-1))/dx;
       drho_dx=(rho(i)-rho(i-1))/dx;
       da_dx=(log(a(i))-log(a(i-1)))/dx;
       dT_dx=(T(i)-T(i-1))/dx;
       
      % continuity equation
       drho_dt_c(i)=-rho(i)*dv_dx-rho(i)*v(i)*da_dx-v(i)*drho_dx;
       
      % momentum equation
       dv_dt_c(i)=-v(i)*dv_dx-1/gamaa*(dT_dx+(T(i)/rho(i))*drho_dx);
      
      % energy equation 
       dT_dt_c(i)=-v(i)*dT_dx-(gamaa-1)*T(i)*(dv_dx+v(i)*da_dx); 
     
       end
   
   % compute the average time derivative:-
    drho_dt_av=0.5*(drho_dt_p+drho_dt_c);
    dv_dt_av=0.5*(dv_dt_p+dv_dt_c);
    dT_dt_av=0.5*(dT_dt_p+dT_dt_c);
    
   % Final solution update:-
     
     for j=2:n-1
      rho(j)=rho_old(j)+drho_dt_av(j)*dt;
      v(j)=v_old(j)+dv_dt_av(j)*dt;
      T(j)=T_old(j)+dT_dt_av(j)*dt;    
     end

   % apply boundary condition:-
   % 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);
     
   % MACH value:-
     M=v./a_conv;
  
   % pressure force value:-
     P=rho.*T;
     
   % Throttle values for parameter:-
   rho_throttle_value(t)=rho(((n-1)/2)+1);
   T_throttle_value(t)=T(((n-1)/2)+1);
   P_throttle_value(t)=P(((n-1)/2)+1);
   M_throttle_value(t)=M(((n-1)/2)+1);
   
  end  
    sim_tim=toc();
  
   % mass flow rate:-
    mass_flow_rate=rho.*v.*a;
    
  text=sprintf(\'NON-CONSERVATIVE FLOW \\n NUMBER OF GRID POINTS USED IN NOZZLE LENGTH=%d,\\n NUMBER OF TIME STEPS=%d\',n,t);
   
  if(t==100) 
   figure(1)
   plot(x,mass_flow_rate,\'b\')
   axis([0 3 0.4 0.7]) 
   xlabel(\'length of nozzle (x/L)\')
   ylabel(\'mass flow rate(\\rho*a*v)/(\\rho_0*a_0*A)\') 
   title(text)
 end
 
   if(t==400) 
   figure(2)
   plot(x,mass_flow_rate,\'b\')
   axis([0 3 0.4 0.7]) 
   xlabel(\'length of nozzle (x/L)\')
   ylabel(\'mass flow rate(\\rho*a*v)/(\\rho_0*a_0*A)\') 
   title(text)
 end
 
   if(t==700) 
   figure(3)
   plot(x,mass_flow_rate,\'b\')
   axis([0 3 0.4 0.7]) 
   xlabel(\'length of nozzle (x/L)\')
   ylabel(\'mass flow rate(\\rho*a*v)/(\\rho_0*a_0*A)\')
   title(text) 
 end
 
   if(t==1400) 
   figure(4)
   plot(x,mass_flow_rate,\'b\')
   axis([0 3 0.4 0.7]) 
   xlabel(\'length of nozzle (x/L)\')
   ylabel(\'mass flow rate(\\rho*a*v)/(\\rho_0*a_0*A)\') 
   title(text)
   end
   
   if(t==1400)
    figure(5)
    subplot(4,1,1)
    plot(rho_throttle_value,\'b\')
    xlabel(\'Number of iteration\')
    ylabel(\'density (\\rho)/(\\rho_0)\')
    title(text)
    
    subplot(4,1,2)
    plot(T_throttle_value,\'r\')
    xlabel(\'Number of iteration\')
    ylabel(\'temperature(T/T_0)\')
    title(text)
    
    subplot(4,1,3)
    plot(P_throttle_value,\'g\')
    xlabel(\'Number of iteration\')
    ylabel(\'pressure (P/P_0)\')
    title(text)
    
    subplot(4,1,4)
    plot(M_throttle_value,\'y\')
    xlabel(\'Number of iteration\')
    ylabel(\'mach number(M)\')
    title(text)
        
    end
   

OUTPUT :-

1.) NON-CONSERVATIVE :-

2.) CONSERVATIVE :-

HERE I OBSERVED THAT THE VALUE OF MASS FLOW IS CONVERGING IN NON-CONSERVATIVE FORM BETTER THAN THE CONSERVATIVE FORM. 

IN BOTH FORM, THE VALUE WAS TOO LARGE FOR LESS TIME STEP i.e UPTO 400 WHICH CAN BE EASILY SEEN IN THE PLOTTED GRAPH.

AT 700 TIME STEP, CONVERGENCE WAS MORE IN CASE OF NON-CONSERVATIVE STATE THAN THE CONSERVATIVE STATE. 

AT 1400 TIME STEP, CONVERGENCE WAS HIGHLY ATTAINED THERE IN CONSERVATIVE WHILE IT WAS  STILL VARYING VERY SMALL MARGINALLY IN CASE OF NON-CONSERVATIVE FORM.

SO IT IS CONCLUDED THAT NON-CONSERVATIVE HAS ADVANTAGE OF ATTAINING CONVERGENCE IN LESS TIME STEP THAN THE CONSERVATIVE FORM.

AND NON-CONSERVATIVE TERM HAVE ONE MORE ADVANTAGE THAT IT IS MORE EXACT NEARER TO THE EXACT VALUE OF MASS FLOW.

WHEN BOTH CASE IS TESTED BY THE EFFECT OF VARYING GRID POINTS, IT IS NOTICED THAT THERE IS HUGE DIFFERENCE DIFFERENCE BETWEEN BOTH CASE WHEN COMPARED AT MINIMUM TIME STEP. WHERE AS, IN CASE OF HIGH TIME STEP OR AT AND AFTER THE TIME STEP VALUE AT WHICH CONVERGENCE OCCUR, IT HAVE VERY LESS DIFFERNCE i.e AROUND ORDER OF `10^-2`.

SO IT IS CONCLUDED THAT THERE IS NO EFFECT OF GRID SIZE VARIATION ONCE THE CONVERGENCE IT IS ACHIEVED. 

HERE WHEN BOTH IS COMPARED IN TERMS OF SIMULATION TIME, IT IS NOTICED THAT NON-CONSERVATIVE CONSUMES LESS TIME FOR SIMULATION WHEN COMPARED TO THE CONSERVATIVE FORM.

AT 31 GRID POINTS:-