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

AIM:-

• To simulate the isentropic flow through a quasi 1D subsonic-supersonic nozzle.
• Derive both the conservation and non-conservation forms of the governing equations and solve them using the MacCormack's technique.
• To determine the steady-state temperature distribution for the flow-field variables and investigate the difference between the two forms of governing equations by comparing their solutions.
• To plot:-

1. Steady-state distribution of primitive variables inside the nozzle.
2. Time-wise variation of the primitive variables at throat for both of the form.

THEORY:-

• Subsonic-Supersonic Isentropic Convergent Divergent Nozzle Flow

Fig.1.0 Subsonic-Supersonic Isentropic Convergent Divergent Nozzle Flow

Source:- COMPUTATIONAL FLUID DYNAMICS, The Basics with Applications, John D. Anderson, Jr. (page.no.285)

The above Fig.1.0 Shows that Subsonic-Supersonic Isentropic Convergent Divergent Nozzle Flow consist of Lenght of L, having two sections Convergent and Divergent seperated by the Throat.

For Fluid flow through nozzle, the fluid flow variables mentioned in the form of Non Dimensional quantity for Numerical values becomes Eaiser to present and Compare.

• Flow field Variables in Non dimensional form

`T'=(T/T^0)`

`rho'=(rho/rho^0)`

`X=(X/L)`

`V'=(V/alpha_0)`

Where, `alpha_0=`Speed of Sound.

`alpha_0=sqrt(gamma*RT_0)`

`t'=(t/(L/alpha_0))`

`A'=(A/A^star)`

Where, `A^star=`throat Area

• Governing Equations

A. Non Conservative Form

• Continuity Equation:-

`(delrho')/(delt')=-rho'*(delV')/(delx')-rho'*V'*((del*ln(A'))/(delx'))-V'*(delrho')/(delx')`

• Mommentum Equation

`(delV')/(delt')=-V'*(delV')/(delx')-(1/gamma)*((delT')/(delx')+(T')/(rho')*(delrho')/(delx'))`

• Energy Equation

`(delT')/(delt')=-V'*(delT')/(delx')-(gamma-1)*T'*((delV')/(delx')+V'*(del*(ln(A')))/(delx'))`

B. Conservative Form

• Continuity Equation:-

`(del(rho'*A'))/((delt'))+(del(rho'*A'*V'))/((delx'))=0`

• Mommentum Equation

`(del(rho'*A'*V'))/((delt'))+(del(rho'*A'*V'^2+1/gamma*rho'*A'))/((delx'))=1/gamma*rho'*(delA')/(delx')`

• Energy Equation

`(del(rho'*((e')/(gamma-1)+gamma/2*V'^2)*A'))/((delt'))+(del(rho'*((e')/(gamma-1)+gamma/2*V'^2)V'*A'+rho'*A'*V'))/(delx')=0`

The above equation can be written as,

`(delU_1)/(delt')=(delF_1)/(delx')`

`(delU_2)/(delt')=(delF_2)/(delx')+J_2``(delU_3)/(delt')=(delF_3)/(delx')`

`F_1=U_2`

`F_2=(U_2^2)/(U_1)+(gamma-1)/(gamma)*(U_3-gamma/2*(U_2^2)/(U_1))`

`F_3=(gamma*U_2*U_3)/(U_1)-(gamma*(gamma-1))/(2)*(U_2^3)/U_1^2`

`J_2=1/gamma*rho'*T'*(delA')/(delx')`

• Meshing

Fig.2.0 Grid point along Nozzle length

Source:- COMPUTATIONAL FLUID DYNAMICS, The Basics with Applications, John D. Anderson, Jr. (page.no.299)

Calculations of Timesteps:-

`Deltat_i=C*(Deltax)/(a_1+V_i)`

`a=sqrtT=`Local Speed of Sound;

C=0.5= Courent Number;

`Deltat=min(Deltat_1,Deltat_2,...Deltat_n)`

Discretization (macormack's method):

1.Predictor method

`(delf)/(delt_p)=`Forward difference Spce term

2. Solution Update

`f^(star)=f_i+(delf)/(delt_p)*Deltat`

3. Corrector Step

`(delf^star)/(delt_c)=`Backward Difference at space term.

`(delf)/(delt_(avg))=0.5*((delf)/(delt_p)+(delf)/(delt_c))`

4.Final Solution Update

`f^2=f_i+(delf)/((delt_(avg)))*(Deltat)`

• Boundary Condition:-

Fig.3.0 Boundary condition Subsonic inflow and Supersonic outflow

Source:- COMPUTATIONAL FLUID DYNAMICS, The Basics with Applications, John D. Anderson, Jr. (page.no.304)

Above fig 3.0 shows that Flow is unsteady and invicid, for 1D unsteady flow their exist 2 real Characteristics line.

1. If characteristics line flow from Left to Right , Drichilet boundary condition Applied Which help to calculate domain of influence.
2. If characteristics line flow from Right to Left Neumannboundary condition Applied Which help to calculate values of boundary by using inner nodes. boundary condition should be Floated.

• Important Formulae:-

1. Nozzle Area(A)=`1+2.2*(x-1.5)^2`

Where, `0lexle3`

2.Pressure (P)=`rho*T`

3.Mach Number =`(V)/(sqrtT)`

4.Mass flow rate=`rho*A*V`

MATLAB CODE:-

MAIN CODE:-

% main code 
clear all
close all
clc

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


tic
[t,rho,p,v,mach_no,mass_flow_rate,t_throat,rho_throat,p_throat,v_throat,mach_throat]=non_conservative_form(n,x,dx);
time_1=toc;
tic
[t_c,rho_c,p_c,v_c,mach_no_c,mass_flow_rate_c,t_throat_c,rho_throat_c,p_throat_c,v_throat_c,mach_throat_c]=conservative_form(n,x,dx);
time_2=toc;

figure(1)
subplot(4,1,1)
plot(x,t,'b','LineWidth',1.3)
ylabel('Temperature')
title('Non conservative form: Distribution of Variables along the nozzle length')
subplot(4,1,2)
plot(x,rho,'b','LineWidth',1.3)
ylabel('Density')
subplot(4,1,3)
plot(x,p,'g','LineWidth',1.3)
ylabel('Pressure')
subplot(4,1,4)
plot(x,mach_no,'r','LineWidth',1.3)
xlabel('Nozzle length')
ylabel('mach-no')


figure(2)
legend('first time step','50th time step','100th time step','200th time step','700th time step','1400th time step')
title('non conservative form: variation of mass flow rate distribution at different time steps')
xlabel('Nozzle Length')
ylabel('mass flow rate')


figure(3)
subplot(4,1,1)
plot(t_throat,'b','LineWidth',1.3)
ylabel('temperature')
title('non conservative form: timewise variations of prmiatevive variables at throat')
subplot(4,1,2)
plot(rho_throat,'b','LineWidth',1.3)
ylabel('Density')
subplot(4,1,3)
plot(p_throat,'g','LineWidth',1.3)
ylabel('Pressure')
subplot(4,1,4)
plot(mach_throat,'r','LineWidth',1.3)
xlabel('Time iterations')
ylabel('mach-no')

figure(4)
subplot(4,1,1)
plot(x,t_c,'b','LineWidth',1.3)
ylabel('temperature')
title('conservative form: Distribution of Variables along the nozzle length')
subplot(4,1,2)
plot(x,rho_c,'b','LineWidth',1.3)
ylabel('Density')
subplot(4,1,3)
plot(x,p_c,'g','LineWidth',1.3)
ylabel('Pressure')
subplot(4,1,4)
plot(x,mach_no_c,'r','LineWidth',1.3)
xlabel('Nozzle Length')
ylabel('mach-no')

figure(5)
legend('first time step','50th time step','100th time step','200th time step','700th time step','1400th time step')
title('conservative form: variation of mass flow rate distribution at different time steps')
xlabel('Nozzle Length')
ylabel('mass flow rate')

figure(6)
subplot(4,1,1)
plot(t_throat_c,'b','LineWidth',1.3)
ylabel('temperature')
title('conservative form: timewise variations of prmiatevive variables at throat')
subplot(4,1,2)
plot(rho_throat_c,'b','LineWidth',1.3)
ylabel('Density')
subplot(4,1,3)
plot(p_throat_c,'g','LineWidth',1.3)
ylabel('Pressure')
subplot(4,1,4)
plot(mach_throat_c,'r','LineWidth',1.3)
xlabel('Time iterations')
ylabel('mach-no')

figure(7)
plot(x,mass_flow_rate,'LineWidth',1.3)
hold on
plot(x,mass_flow_rate_c,'m','LineWidth',1.3)
hold on
legend('Non conservative_form','Conservative form')
xlabel('Nozzle length')
ylabel('mass flow rate')
title('comparison of Normalized mass flow rate distributions of both forms')

FUNCTION CODE:-

For Non-conservative Form:-

function [t,rho,p,v,mach_no,mass_flow_rate,t_throat,rho_throat,p_throat,v_throat,mach_throat]=non_conservative_form(n,x,dx)

% initial profiles
rho=1-0.3146*x;
t=1-0.2314*x;
v=(0.1+1.09*x).*(t).^0.5;

gamma=1.4;

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

% time steps
nt=1400;

dt=min(0.5*dx./(v+sqrt(t)));

% outer time step

for k=1:nt
    % copy of variables 
    rho_old=rho;
    v_old=v;
    t_old=t;

    % predictor method
    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;

        % continuity equation
        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 method
    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;

        % continuity equation
        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
        % average 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 upodate
        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

        % 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);

        % calculation pressure, mach no and mass flow rate

        p=rho.*t;
        mach_no=v./sqrt(t);
        mass_flow_rate=rho.*a.*v;

        % calculation pressure, mach no and mass flow rate at throat

        p_throat(k)=p(16);
        v_throat(k)=v(16);
        mach_throat(k)=mach_no(16);
        t_throat(k)=t(16);
        rho_throat(k)=rho(16);

        % variation of mass flow rate of distribution inside  the nozzle at
        % different time steps
        figure(2)
        if k==1
            plot(x,mass_flow_rate,'LineWidth',1.3)
            hold on
        elseif k==50
                plot(x,mass_flow_rate,'r','LineWidth',1.3)
                hold on
        elseif k==100
                plot(x,mass_flow_rate,'b','LineWidth',1.3)
                hold on
        elseif k==200
                plot(x,mass_flow_rate,'cyan','LineWidth',1.3)
                hold on
        elseif k==700
                plot(x,mass_flow_rate,'m','LineWidth',1.3)
                hold on
        elseif k==1400
                plot(x,mass_flow_rate,'g','LineWidth',1.3)
                hold on
        end
end
end
      
       

For Conservative Form:-

%%Function file 

function [t_c,rho_c,p_c,v_c,mach_no_c,mass_flow_rate_c,t_throat_c,rho_throat_c,p_throat_c,v_throat_c,mach_throat_c]=conservative_form(n,x,dx)

% initial profile

for i=1:n
    if (0<=x(i) && 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)
        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;
gamma=1.4;
v_c=0.59./(rho_c.*a);

% solution vector
u1=rho_c.*a;
u2=rho_c.*a.*v_c;
u3= rho_c.*a.*((t_c./(gamma-1))+(gamma/2).*v_c.^2);


% time steps
nt=1400;
dt=min(0.5*dx./(v_c+sqrt(t_c)));

% Outer time steps
for k=1:nt

    % copy of u1,u2,u3
    u1_old=u1;
    u2_old=u2;
    u3_old=u3;

    % Flux vectors
    
    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));

   % predictor method
   for j=2:n-1
       j2 = (1/gamma)*rho_c(j)*t_c(j)*((a(j+1)-a(j))/dx);
       % continuity equation
       du1dt_p(j) = -(f1(j+1)-f1(j))./dx ;
       % mommentum equation
       du2dt_p(j) = -(f2(j+1)-f2(j))./dx + j2 ;

       % Energy equation
        du3dt_p(j) = -(f3(j+1)-f3(j))./dx ;

     % solution update to be done
       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 primitive variables
   rho_c=u1./a;
   v_c=u2./u1;
   t_c = (gamma-1).*((u3./u1) - (gamma./2).*(v_c).^2);

   % updating flux vectors
    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));


   % corrector Method
   for j=2:n-1

      j2 = (1/gamma)*rho_c(j)*t_c(j)*((a(j)-a(j-1))/dx);
      % Correctors Step derivative
      du1dt_c(j) = -(f1(j)-f1(j-1))./dx ;
      du2dt_c(j) = -(f2(j)-f2(j-1))./dx + j2 ;
      du3dt_c(j) = -(f3(j)-f3(j-1))./dx ;
   end

       % average 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 m=2:n-1
       u1(m)=u1_old(m)+du1dt(m)*dt;
       u2(m)=u2_old(m)+du2dt(m)*dt;
       u3(m)=u3_old(m)+du3dt(m)*dt;
  end

       % Boundary condition
       % Inlet

    u1(1) = rho_c(1).*a(1);
    u2(1)= 2*u2(2)-u2(3);
    u3(1)= u1(1)*((t_c(1)/(gamma-1))+(gamma/2)*v_c(1).^2);

    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 update on primitive variables
       rho_c=u1./a;
       v_c=(u2)./(u1);
       t_c = (gamma-1)*((u3./u1)-(gamma/2)*v_c.^2);

       % calculation of pressure, mach no, and mass flow rate
       p_c=rho_c.*t_c;
       mach_no_c=v_c./sqrt(t_c);
       mass_flow_rate_c=rho_c.*a.*v_c;
       % calculation of pressure, mach no, and mass flow rate at throat
       p_throat_c(k)=p_c(16);
       v_throat_c(k)=v_c(16);
       mach_throat_c(k)=mach_no_c(16);

       t_throat_c(k)=t_c(16);
       rho_throat_c(k)=rho_c(16);

        % variation of mass flow rate of distribution inside  the nozzle at
        % different time steps

        figure(5)
        if k==1
            plot(x,mass_flow_rate_c,'LineWidth',1.3)
            hold on
        elseif k==50
                plot(x,mass_flow_rate_c,'r','LineWidth',1.3)
                hold on
        elseif k==100
                plot(x,mass_flow_rate_c,'b','LineWidth',1.3)
                hold on
        elseif k==200
                plot(x,mass_flow_rate_c,'cyan','LineWidth',1.3)
                hold on
        elseif k==700
                plot(x,mass_flow_rate_c,'m','LineWidth',1.3)
                hold on
        elseif k==1400
                plot(x,mass_flow_rate_c,'g','LineWidth',1.3)
                hold on
        end
end
end

STEP-WISE EXPLAINATION:-

Step-1:- First of all Initialize the value of Grid Number(n), x and dx values.

Step-2:- Then Assign the values of Variables(initial profiles) for Conservative and Non Conservative form by Creating Function Seperatly.

Step-3:- For mac-cormack method Outer time step loop is created and it was initiated.

Step-4:- Predictor method and Corretor method carried out and also final solution update also done for both methods.

Step-5:- Boundary conditions are applied for Inlet and Outlet, the calculate the Pressure, mach number and mass flow rate.

Step-6:- the values of primitive are assigned for calculating variables at each timestep at throat, then, mass flow rate are plotted at different different time steps.

Step-7:- then same proccess carrid out for Conservation form except Flux vectors and Solution Vectors (F1,F2 andF3) and (U1,U2 and U3) Respectively.

Step-8:- In mac-cormack method Flux vectors and Solution Vectors used and then primitive variables are updated.

Step-9:- above procedure repeated up to 1400 time-steps and Function return value in main code.

Step-10:- Then Respective plots are Plotted by using figure and subplot Commands.

OUTPUT:-

NON CONSRTVATIVE FORM:-

• Distribution of Variables (Temperature, Density, Pressure and Mach Number) along the nozzle length.

Fig.1 Distribution of Variables (Temperature, Density, Pressure and Mach Number) along the nozzle length.

Fig.2.Variation of mass flow rate distribution at different time steps

Fig.3.Timewise variations of prmiatevive variables (Temperature, Density, Pressure and Mach Number) at throat.

EXPLAINATION:-

1. In Fig 1. The temperature Density and pressure values decreased as Nozzle lenght increases, but mach Number incresed as Nozzle length Increses.

2. In Fig 2. As mass flow rate is plotted, the mass flow rate does not change at 700th time step to 1400 time step, so we can say almost reached steady state at 700th timestep.

3. In Fig 3. Fluctuations observed at initial time steps, then plot becomes linear which means that reached steady state.

CONSERVATIVE FORM:-

• Distribution of Variables (Temperature, Density, Pressure and Mach Number) along the nozzle length.

Fig.4. Distribution of Variables (Temperature, Density, Pressure and Mach Number) along the nozzle length

• Variation of mass flow rate distribution at different time steps:-

Fig.5. Variation of mass flow rate distribution at different time steps.

• Timewise variations of prmiatevive variables (Temperature, Density, Pressure and Mach Number) at throat:-

Fig.6.Timewise variations of prmiatevive variables (Temperature, Density, Pressure and Mach Number) at throat.

EXPLAINATION:-

1. In Fig 4. The temperature Density and pressure values decreased as Nozzle lenght increases, but mach Number incresed as Nozzle length Increses.

2. In Fig 5. As mass flow rate is plotted, the mass flow rate does not change at 700th time step to 1400 time step, so we can say almost reached steady state at 700th timestep.

3. In Fig 6. Fluctuations observed at initial time steps, then plot becomes linear which means that reached steady state.

Comparison of Mass Flow rate:-

Fig.7. Comparison of Mass Flow rate

EXPLAINATION:

In fig.7. We can see that Better mass flow rate at Conservative form found than Non-conservative form more fluctuation found in Non Conservative Form.

GRID INDEPENDENCE TEST:-

When Grid Size is having n=31 primitive values at throat region of the nozzle as shown in Below table. also checked primitive values for Grid Size 61, in that case we found no any change in values, So the values for Conservative and Non-conservative are Almost closer. so we found Change in grid size does not affect the values of primitive varibles.

Grid Independence test at nozzle throat:-

	No. of grid points	Temerature	Pressure	Density	mach Number	Simulation Time
Non Conservative	31	0.8365	0.5342	0.6387	0.9994	1.2915
Conservative	31	0.8411	0.5487	0.6524	0.9780	0.9224
Non Conservative	61	0.8354	0.5327	0.6377	0.9998	0.9662
Conservative	61	0.8357	0.5340	0.6390	0.9978	0.9269

Table.1. Grid Independence test at nozzle throat.

Observations:-

1. Simulation time is lesser in Conservation Form than Non- conservation in both cases of Number of Grids.

2. No any Difference in Results of Conservative and Non-Conservative form as primitive varibles values concerns.

CONCLUSION:-

• In MATLAB The Isentropic Flow through 1D quasi subsonic -Supersonic Nozzle Simulated.
• The Governing Equations are solved for Conservation And Non-Conservation Form by using MacCormack Method.
• Conservation and Non-conservation Form found Grid Indepened.
• Mass flow rate Distribution are plotted and We found better mass flow rate In Conservation Form.

REFERENCES:-

BOOK:- COMPUTATIONAL FLUID DYNAMICS, The Basics with Applications, John D. Anderson, Jr.