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

ANALYSIS OF DE-LAVAL NOZZEL

The solution for a 1D supersonic nozzle can be obtained by MACCORMACK method. The steps involved in analysis of such a problem are described below

• Physical description of the fluid flow
• CFD intermediate steps
• CFD final solution

Step 1 PHYSICAL DESCRIPTION OF FLUID FLOW

we have a C-D nozzle given as a problem statement let us understand how the cd nozzle accelerates the flow from subsonic to supersonic and the conditions required for them .

figure 1: c-d nozzle

the convergent divergent nozzle accelerates the flow from stagnation condition to supersonic condition.

In this project we are calculating for two forms of equations

• Conservative form
• Non conservative form

in this section lets first discuss about the non conservative form as it simply uses the primitives for analysis .

equation 1: general equation of continuity , momentum and energy equations .

these are the set of equation that has to solved for all the grids before that lets non dimentionalise them to get a set of equtaions for continuity momentum and energy .

equation 2: the primitive form or non dimentional form of the above mentioned three equation.

these are the primmitive non conservation equations .we will be solving them numerically how to arrive to this method is fairly simple we just use the non dimentional gas dynamic properties.

Here we employ primitives that is all the general form of continuity ,momentum and energy equations are non dimentionalized as done in Anderson . the main aim of this project is to examine the difference in the two forms when it comes to CFD solutions

Step 2 CFD INTERMEDIATE STEP

The physical flow properties and equations pertaining to it are derived in their non dimentional form.

Maccromack method

we will set up a time-marching fimted1fference solut10n for the quasi-one-dimensional nozzle flow.we have the problem defined as a function of x direction alone since it is quasi one dimentional flow we solve for one dimentional flow field with time marcing solution .the set up is the same foe both conservative and non conservative form just the equations change a little .

• predictor step we set up the spatial derivatives as forward differences. Also, to reduce the complexity of the notation, we will drop the use of the prime to denote a dimensionless variable. all variables are the nondimensional variables.

equation 3 : the time marching solving for the predictor step at time (t).

Once these three equations are calculated we need to update them for a small increment of time at (t+dt) as the above calculated values are for the time step (t)

equation 4 : the solution update for the primitives at time (t+dt).

once we have the predicted values we update the parameters of study which are later used for the corrected step ,there is no differnece in the merthod for predictor and corrector step only thing is in corrector step we use updated values and we use the average value od both predictor and corrector values and calculate the incremental time step of parameters one of the important point to be noted is in corrected step wwe use backward difference method for space terms decribed in detail in computaional fluid dynamics anderson textbook.the corrected value for (t+dt)time step is given by the below formaulae .

equation 5 : the final solution update for the primitives using the average values at (t+dt).

we can see that we use the average values of the properties calculated from predictor and corrected step which are updated for (t+dt)step once this step is completed then the solution is fairly arrived we need to calculate the values at every grid points along the x axis.

These are the updated values after every time step and corrected values .this is the maccromack method of solving the problem 1d quasi supersonic nozzle .

lets move on to an interesting part off calculating the time step we can clearly see that we have an hyperbolic equation for the parametes and we are using an explicit method here we should be carefull in selecting the time step which is given by

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

where C=coutrant number

dx=space step

a= local speed of sound

v= local velocity

these values are calculated at each and every grid points and the minimum value along all the grid point is taken and then the minimum of that time step is used for calculating IN OUR ANALYSIS WE ARE USING A GUESS VALUE OF dt TO UNDERSTAND THE EFFECT OF CHANGE IN dt .

we arrive at dt=0.002 the solution is stable therfore its better to change the value and understand the effect on stability

BOUNDARY CONDITIONS

DEFINING BOUNDARY CONDITIONS

Since the geometry is 1D we discretize along x axis with N grid points ,the flow enters at the resorvior i=1 and it will be the intake and n will be the outflow of the nozzle

We fix the density and temperature as unity and allow velocity to float as one of the parameter is allowed to float to maintain mass flow rate as it is allowed to adjust itself .

We obtain velocity by using the slope function and v1=2*(v2-v3);using the extrapolation method

figure 2: velocity vs length (dx) for linaer extrapolation

This is the linear extrapolation curve .Linear extrapolation method figure.

Supersonic outflow condition here we are allowing all the parameters to float we cannot constrain any parameter of nozzle hence we keep all the parameters to float.

equation 7: the outflow boundary conditions

N is the last node in the grid .

INITIAL CONDITIONS

The most important is the design of the nozzle we have to define the area parameter that is non dimentional area ratio A/A*;

A/A* = I + 2.2*(x - I.5)^2;

X varies from 0 to 3 units . The closer the initial conditions are to the final steady-state answer, the faster the time-marching procedure will converge, and hence the shorter will be the computer execution time hence we have a choose a proper initial conditions for the problem statement from the subsonic region to the supersonic regime since the pressure and density decreses along the x axis we have choose a value such a function that these value decrease with space coordinates.

equation 8:the inital condition for the nozzle.

LETS MAKE IT MORE INTERESTING JUST CHANGE THE INITIAL PROFILE AND OBSERVE THE CHANGES THE SOLUTION REMAINS THE SAME UNTILL THE VALUES OF DENSITY AND TEMPERATURE DECREASES OVER THE LENGTH AND VELOCITY INCREASES OVER LENGTH IF THE CONDITION DOES NOT SATISFY THE SOLUTION DOES NOT CONVERGE .

***************************main program *******************************************

clear all
close all 
clc
 n=31;
 x=linspace(0,3,n);
 dx=3/(n-1);
 nt = 1400;
 throat = ((n-1)/2)+1;
 % calling the function 
 %conservative form 
 [zc,tempc,ms,mc,vc,pc,time1,z_t,tempc_t,pc_t,mc_t,vc_t,ms_t] = nonconservativeformnope(n,x,dx,nt,throat);
 %nonconservative

% plot(x,pc)
% plot(x,tempc)
% plot(x,m)
%plot(x,msc)
[z,v,t,msc,m,p,time2,zz_t,temp_t,p_t,m_t,v_t,msc_t]=conservativeform(n,x,dx,nt,throat);
exactp=[1 0.997 0.996 0.992 0.988 0.982 0.974 0.963 0.948 0.926 0.895 0.853 0.795 0.722 0.633 0.534 0.434 0.343 0.264 0.201 0.153 0.117 0.090 0.070 0.055 0.044 0.035 0.029 0.024 0.020 0.014]

%for 1400time steps lets find the difference 
diff_p=abs(pc-exactp)
figure(1);
subplot(4,1,1)
plot(x,pc,'k','LineWidth',1);
ylabel('pressure');
title('parameters for nonConservative equations');
grid on;
subplot(4,1,2)
plot(x,zc,'k','LineWidth',1);
ylabel('Density');
grid on;
subplot(4,1,3)
plot(x,tempc,'k','LineWidth',1);
ylabel('Temperature');
grid on;
subplot(4,1,4)
plot(x,mc,'k','LineWidth',1);
ylabel('mach no');
xlabel('dimonstionless Distance of nozzle');
grid on;

figure (2)
plot(x,pc,'*',x,exactp,'r','LineWidth',1.2)
title ('exactvsnumerical(non conservative form)')


legend('numerical','exact')
xlabel('NON DIMENTIONAL DISTANCE')
ylabel('presuure values')

figure (3)
yyaxis right
plot(x,mc,'b','LineWidth',1)
title('density and mach no variation')
hold on 
yyaxis left
plot(x,zc,'g','LineWidth',1)
legend('density','mach no')
xlabel('NON DIMENTIONAL DISTANCE')
grid on


% plotting for conservative 
figure(4);
subplot(4,1,1)
plot(x,p,'k','LineWidth',1);
ylabel('pressure');
title('parameters for Conservative equations');
grid on;
subplot(4,1,2)
plot(x,z,'k','LineWidth',1);
ylabel('Density');
grid on;
subplot(4,1,3)
plot(x,t,'k','LineWidth',1);
ylabel('Temperature');
grid on;
subplot(4,1,4)
plot(x,m,'k','LineWidth',1);
ylabel('mach no');
xlabel('dimonstionless Distance of nozzle');
grid on;

figure(5)
plot(x,msc,x,ms,'r','LineWidth',1.2)
title('massflow')
xlabel('dimonstionless Distance of nozzle');
legend('mass flow conservative form','mass flow non conservative form')
grid on 

figure (6)
subplot(4,1,1)
plot(x,pc,'*',x,p,'r','LineWidth',1.2)
xlabel('dimonstionless Distance of nozzle');
ylabel('pressure value');
subplot(4,1,2)
plot(x,tempc,'*',x,t,'r','LineWidth',1.2)
xlabel('dimonstionless Distance of nozzle');
ylabel('temperature value');
subplot(4,1,3)
plot(x,zc,'*',x,z,'r','LineWidth',1.2)
xlabel('dimonstionless Distance of nozzle');
ylabel('density value');
hold on
subplot(4,1,4)
plot(x,mc,'*',x,m,'r','LineWidth',1.2)
xlabel('dimonstionless Distance of nozzle');
ylabel('mach value');
legend(' conservative form',' non conservative form')
title('comparitive study of conservative and non conservative form  ')
grid on 

figure(7)
plot(m,p)
hold on
plot(m,z)
hold on
plot(m,t)
grid on
xlabel('mach no');
ylabel('parameter value');
legend(' pressure','density ','temperature')
title('variation of properties with mach number')


figure(8)
subplot(4,1,1);
plot(1:nt,zz_t);
hold on
plot(1:nt,z_t, 'color', 'r');
grid on
ylabel('Density');
legend('Conservative form' , 'Non-conservativeform');
title('Study of density variation at each time step at throat');

subplot(4,1,2);
plot(1:nt,temp_t);
hold on
plot(1:nt,tempc_t);
grid on
ylabel('Temperature');
legend('Conservative form' , 'Non-conservativeform');
title('Study of tempertaure variation at each time step at throat');

subplot(4,1,3);
plot(1:nt,p_t);
hold on 
plot(1:nt,pc_t);
grid on
ylabel('Pressure');
legend('Conservative form' , 'Non-conservativeform');
title('Study of pressure variation at each time step at throat');


subplot(4,1,4);
plot(1:nt,m_t);
hold on
plot(1:nt,mc_t);
grid on
xlabel('Time step')
ylabel('Mach Number');
legend('Conservative form' , 'Non-conservativeform');
title('Study of Mach number variation at each time step at throat');

in the main program we are calling two different functions that solve for conservative and non conservative forms of equations and plot the paramters of interest we give some arguments for both the functions .and the functions incorporates macromac method for solving which is detailed above

now lets call one function nonconservative form and plot the results and compare the plots

function [z,temp,ms,m,v,p,time1,z_t,tempc_t,pc_t,mc_t,vc_t,ms_t]=nonconservativeformnope(n,x,dx,nt,throat)
%x=linspace(0,3,31)
%defining the area of the nozzle

a=1+(2.2*(x-1.5).^2);
plot(a);
%defining paparemts

z=1-(0.3146*x);
temp=1-(0.2314*x);
v=(0.1+(1.09*x)).*temp.^0.5;

cfl=0.5;
gamma =1.4;
%pareameter of studty
dt = 0.02;

%main program 
tic 
for k=1:nt
    
    %we are taking a copy before execution
    zold=z;
    vold=v;
    tempold=temp;
    
    %STEP 1 predictor step 
    for j=2:(n-1)
        
        %discritization of equations 
        dzdx=(z(j+1)-z(j))/dx;
        dvdx=(v(j+1)-v(j))/dx;
        dtdx=(temp(j+1)-temp(j))/dx;
        dadx=(log(a(j+1))-log(a(j)))/dx;
        
    %equations to solve
    dzdt_p(j)= -(z(j)*dvdx)-z(j)*v(j)*dadx-v(j)*dzdx;
    
    dvdt_p(j)= -v(j)*dvdx-(1/gamma)*(dtdx+(temp(j)/z(j))*dzdx);
    
    dtdt_p(j)= -v(j)*dtdx-(gamma-1)*temp(j)*(dvdx+v(j)*dadx);
    
    %update step
    z(j)=z(j)+dzdt_p(j)*dt;
    v(j)=v(j)+dvdt_p(j)*dt;
    temp(j)=temp(j)+dtdt_p(j)*dt;
    end
    %STEP 2 corrector step
    %BACK WARD DIFFERNECE
     for j=2:(n-1)
        
        %discritization of equations 
        dzdx=(z(j)-z(j-1))/dx;
        dvdx=(v(j)-v(j-1))/dx;
        dtdx=(temp(j)-temp(j-1))/dx;
        dadx=(log(a(j))-log(a(j-1)))/dx;
        
    %equations to solve
    dzdt_c(j)= -(z(j)*dvdx)-z(j)*v(j)*dadx-v(j)*dzdx;
    
    dvdt_c(j)= -v(j)*dvdx-(1/gamma)*(dtdx+(temp(j)/z(j))*dzdx);
    
    dtdt_c(j)= -v(j)*dtdx-(gamma-1)*temp(j)*(dvdx+v(j)*dadx);
     end
    %update step
   dzdt=0.5*(dzdt_p+dzdt_c);
   dvdt=0.5*(dvdt_p+dvdt_c);
   dtdt=0.5*(dtdt_p+dtdt_c);
   
   %final values
   for j=2:(n-1)
       z(j)=zold(j)+dzdt(j)*dt;
       v(j)=vold(j)+dvdt(j)*dt;
       temp(j)=tempold(j)+dtdt(j)*dt;
   end
   
 
    %Boundary conditions
        
       v(1) = 2*v(2) - v(3);
       z(n) = 2*z(n-1) - z(n-2);
       v(n) = 2*v(n-1) - v(n-2);
       temp(n) = 2*temp(n-1) - temp(n-2);
       
       %parameters of calculation 
        p=z.*temp;
        
        m= v./(temp.^0.5)
      ms=z.*a.*v;
      %throat values 
        z_t(k) = z(throat);
    tempc_t(k) = temp(throat); 
    pc_t(k) = p(throat);
    mc_t(k) = m(throat);
    vc_t(k) = v(throat);
    ms_t(k) = ms(throat)
    
    figure(9)
    if k == 50
        plot(x,ms)
        hold on
        grid on
    elseif k==100
        plot(x,ms)
    elseif k==150
        plot(x,ms)
    elseif k==200
        plot(x,ms)
     elseif k==700
        plot(x,ms)
    
 xlabel('Length of Nozzle ');
ylabel('Mass flow rate ');
title('Non dimensional mass flow rate at different Time steps for non conservative form');
legend('50 Time steps','100 Time steps','150 Time steps','200 Time steps','700 Time step');
      
      
     end
     toc
     time1=toc
     
     
end
   
        
     
    

PLOTS AND DISCUSSIONS

under this section we are understanding the plots and the parametrs variation in the nozzle .

PLOT 1

in this plot we are considering four parametrs that affect the function of the nozzle pressure density and temperature and mach number here we are not considering velocity explicitly as they are incorporated into mach number .

the plot is exactly the required plot we get a smooth curve as the function we are studying is hyperbolic form we can see the function variation in the same manner along the length of the nozzle if any one of this parameter goes wrong the output that we get for the other parameters completely change ,it mainly depends on the initial condition we feed to the function for density temperature and velocity ,the pressure value is just the output as we use the ideal gas state equation .as mentioned we can play with the initial conditions.

PLOT 2

This plot is just to understand the difference in the value calculated numerically and analalytical solution and mainly we are taking the values of pressure and the length of the nozzle u may wonder why pressure their is a reason for it pressure is the product of density * temperature so we can see that error will multipy hence the difference will be more .apart from this we can see the difference in pressure is of the order (1e-5) therefore the residuals are very less for non conservative form .

PLOT 3

this is just a verification plot to observe how the parameters vary with mach number so that we can conclude if the solution is moving towards stability with increase in mach number we should get the plot such that the pressure density and temperature decrease with mach number .

PLOT 4

the figure on the left is from anderson (computational fluid dynamics) and one on the right is the plot we get from the non conservative equation solution .This comparition is only for understanding the flow physics and comparitive study .

*************************************************************************************************************************

THE NEXT FORM OF EQUATION IS THE CONSERVATIVE FORM

equation 9: the gernal form of naviour stokes equation

this is the general form of continuity equation lets first change this form into the primitive form anyways we wont use the primitive form ,lets see the procedure .

equation 10: the primitives having the non dimontionless form .

we have the three equations continuity ,momentum and energy for conservative form . but we are not solving for the primitives we are changing them into solution ,flux vectors and from them we arrive at the values for the primitives .

equation 11: the solution vectors and flux defined in term of primitives .

using these we can change the primary equations and get a set of partial differntials using flux and solution variables

equation 12: the final PDE to solve numerically.

these three will be the equations that we will be solving for conservative form .their is a reson for why we use the solvution and flux variables to solve .we will be solving for the dependent variables rather than the primitives we will be getting the solution for u1,u2 ,u3 directly after the first time step and we have to decode the vales for the primitives from the solution vectors .

equation 13 :decoding the primes from the solution vectors.

here we get the values for density velocity temperature and pressure values after decoding the solution vectors(u1,u2,u3).

the method for solving using the maccromac method is the same their is no change in the method for solving the partial differential equation from predictor to corrector and update steps they are all the same .

so lets move on directly to the program that solves the conservative equations .

function [z,v,t,msc,m,p,time2,zz_t,temp_t,p_t,m_t,v_t,msc_t]=conservativeform(n,x,dx,nt,throat)

 gamma=1.4;
 
 %define geometry
 a=1+(2.2*(x-1.5).^2);
 
 %initial condition for conservative equation from anderson 
 for i=1:n
     if(x(i)>=0 && x(i)<=0.5)
         z(i)=1;
         t(i)=1;
     elseif (x(i)>=0.5 && x(i)<=1.5)
         z(i)=1-0.366*(x(i)-0.5);
         t(i)=1-0.166*(x(i)-0.5);
         elseif (x(i)>=1.5 && x(i)<=3)
         z(i)=0.634-0.3876*(x(i)-1.5);
         t(i)=0.833-0.35166*(x(i)-1.5);
     end 
 end 
 % defining velocity
 v=(0.59)./(z.*a);
 c=0.5;
 dt=0.02;
 
 %since it is conservative form lets start with solution vectors
 u1=z.*a;
 u2=z.*a.*v;
 u3 = u1.*((t./(gamma - 1)) + ((gamma/2).*(v.^2)));
tic
for k=1:nt
    u1old=u1;
    u2old=u2;
    u3old=u3;
    
   %flux terms in solution vectores rather than primitive form 
   f1=u2;
   f2 = ((u2.^2)./u1) + (((gamma - 1)/gamma)*(u3 - ((gamma/2)*((u2.^2)./u1))));
   f3 = ((gamma*u2.*u3)./u1) - (gamma*(gamma-1)*(u2.^3))./(2*(u1.^2));
   
   % macromac method 
   %step 1 predictor method 
   
   for j=2:n-1
       dadx=(a(j+1)-a(j))/dx;
       
       %source term 
       j2=(z(j)*t(j)*dadx)/gamma;
       
       %solving the pure form of conservative equation
       %continuity eqn
       du1dt_p(j)=-(f1(j+1)-f1(j))/dx;
       
       % momentum eqn 
       du2dt_p(j)= -(f2(j+1)-f2(j))/dx + j2;
       
       %energy equation 
       du3dt_p(j)= -(f3(j+1)-f3(j))/dx;
       
       %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
   
   % calculating primitive terms 
   v_p=u2./u1;
   z_p=u1./a;
   t_p= (gamma - 1)*((u3./u1) - ((gamma/2)*(v_p.^2)));
    f1 = u2;
    f2 = ((u2.^2)./u1) + (((gamma - 1)/gamma)*(u3 - ((gamma/2)*((u2.^2)./u1))));
     f3 = ((gamma*u2.*u3)./u1) - (gamma*(gamma-1)*(u2.^3))./(2*(u1.^2));
     
     % step 2 corrector step 
     
    for j=2:n-1
       dadx=(a(j)-a(j-1))/dx;
       
       %source term 
       j2=(z(j)*t(j)*dadx)/gamma;
       
       %solving the pure form of conservative equation
       %continuity eqn
       du1dt_c(j)=-(f1(j)-f1(j-1))/dx;
       
       % momentum eqn 
       du2dt_c(j)= -(f2(j)-f2(j-1))/dx + j2;
       
       %energy equation 
       du3dt_c(j)= -(f3(j)-f3(j-1))/dx;
    end        
        %Average tsolution update
        du1dt = 0.5*(du1dt_p + du1dt_c);
        du2dt = 0.5*(du2dt_p + du2dt_c);
        du3dt = 0.5*(du3dt_p + du3dt_c);
        
        % final values 
        for j = 2:n-1
            u1(j) = u1old(j) + du1dt(j) * dt;
            u2(j) = u2old(j) + du2dt(j) * dt;
            u3(j) = u3old(j) + du3dt(j) * dt;
        end
         %Boundary Conditions
       % same way as non conservative method 
        u1(1) = z_p(1)*a(1);
        u2(1) = 2*u2(2) - u2(3);
        u3(1) = u1(1)*((t_p(1)/(gamma-1)) + ((gamma/2)*v_p(1)^2));
        v_p(1) = u2(1)/u1(1);
        %outlet 
        %extapolation method 
         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 output values 
        z = u1./a;
        v = u2./u1;
        t = (gamma - 1)*((u3./u1) - ((gamma/2)*(v.^2)));
        msc = z.*v.*a;
        m = v./((t).^0.5);
        p = z.*t;
          zz_t(k) = z(throat);
    temp_t(k) = t(throat); 
    p_t(k) = p(throat);
    m_t(k) = m(throat);
    v_t(k) = v(throat);
    msc_t(k) = msc(throat)
    
      figure(10)
    if k == 50
        plot(x,msc)
        hold on
        grid on
    elseif k==100
        plot(x,msc)
    elseif k==150
        plot(x,msc)
    elseif k==200
        plot(x,msc)
     elseif k==700
        plot(x,msc)
         xlabel('Length of Nozzle ');
ylabel('Mass flow rate ');
title('Non dimensional mass flow rate at different Time steps for conservative form');
legend('50 Time steps','100 Time steps','150 Time steps','200 Time steps','700 Time step');

      
         end
    time2 = toc;
end

       

this is the function for the main program we are solving for conservative form .

lets observe the plots for the conservative form alone.

PLOT 1

the above plot consists parameters vs the length of the nozzle for conservative form of equations we have variation of pressure ,density ,temperature and mach number respectively we can see a similar trend compared to the non conservation form .

for better understanding both forms lets plot few results that explains the different forms.

RESULTS AND DISCUSSIONS

1. Steady-state distribution of primitive variables inside the nozzle

we have seen both the plots they have the similar trend we cannot conclude the solution so lets plot them in the same figure .

here we can see clearly that there is no difference in both the forms of equations let it be conservative or non conservative form for the values of primitives they have similar trend and value of difference is less very negligle hence both fall on the same curve . pressure ,density ,temperature and mach number for conservative and non conservative form.

PLOT 2

2. Time-wise variation of the primitive variables

for the time steps of 1400 we are ploting the variation of primitives and their variation over time .

we are just selecting a node and finding the behaviour of the properties we are having a eulers approach toward finding the variation of the properties .

when we have density, temperature ,pressure and mach number at the throat ,which will be the 16th node along the length of the nozzle(x)lets observe the variation of one of the parameter of study lets take pressure and if we see we have two lines one in red the non conservative part and the other in blue the conservative part .when we see the red line we can see for the first 400 time steps we can see a sine function type of variation not exactly sine but the amplitude decreases a damped function and stabilizes after 400 time steps but when we consider the conservative form the fluctuations it undergoes is very less it looks like a constant line through out .

PLOT 3

3. Variation of Mass flow rate distribution inside the nozzle at different time steps during the time-marching process.

both the plots have the mass flow rate at different time steps here we are plotting like that of the anderson(refference text) for compartative study we plot for five different time steps 50,100,150,200,700 in both conservative and non conservative forms of the equations in all the time spets we can see alot of variation along the length as the solution to trying to reach stability respectively lets study the green line thats for the 700th time step we can see that the mass flow has reached a steady state we can see the change over time is negliglible and the slope is zero .so whats basically happening is the equations are solved over time to arrive at a steady stae solution for both conservative and non conservative forms .

PLOT 4

4. Comparison of Normalized mass flow rate distributions of both forms

We can see that from the plot there is no difference between the two forms of equation conservative and non conservative they yeild the same results but when we come to solving using numerical method we need to choose an appropriate method for solving or we may end up at the wrong solution platform and at the same time when we see the plot of mass flow plot we can see that there is variation at inflow and outflow boundaries for non conservative form while the conservative form is almost constant and is approximatly equal to the exact value.

Therfore we can conclude that conservative form of solving is an advantage the reason being this form of equation ensures the mass remains constant because the mass flow is a dependent variable.the results obtains from conserrvative equation yeilds mass flow directly.

In non conservative part the equation give solution for the primitive values and the mass flow is secondary result from the equations . hence we can see the primitives of both the form have same value and the error among them is comparitively negligible and they are combindly ploted as shown in the figure.

At the throat we can see that the non conservative form are more close towards the solution compared to conservative form so we cannot clearly say which is better ,they have their own advantage .

• conservative form is better in conserving mass and mass distribution through the flow field .
• non conservative form leads to smaller residuals which is preffered for numerical analysis hence they do the job better .non conservatibe for 1E-6
• The non conservation form seems to produce slightly more accurate results for the primitive variables, and the conservation form seems to produce slightly more accurate results for the flux variables.

**************************************************************************************************************************************

GRID INDEPENDENCE STUDY

table 1 : grid independence study

here we are changing the number of grid points to understand how it affects the values of the parametrs of our study we know that with increase in grid points the error pertaining to the solution decreases.which is more important to attain ,once we have fixed our geometry and physics we increase the grid size to attain a solution that remains same in all aspects .

***************************************

note -the above code pasted is for 31 grid points and for 61 grid we have change the value of n to 61 and one plot thats for plotting the values for exact and numberical in line 44 of the amin program has to be changed .