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

Objective:

1) To write codes for conservative and non-conservative form of governing equations of 1D supersonic nozzle flow using MacCormack method.

2) Compare the normalized mass flow rate between the conservative forms and non-conservative forms.

3) To find out which method is faster.

THEORY:

We consider the steady, isentropic flow through a convergent-divergent nozzle. The flow at the inlet to the nozzle comes from a reservoir where the pressure and temperature are denoted by p0, and To. Since the cross-sectional area of the reservoir is large so velocity at the initial position is considered to be 0, Thus, Po and T0 are the stagnation values, or total pressure and total temperature respectively.

The flow expands isentropically to supersonic speed at the exit of the nozzle. where the exit pressure, temperature, velocity, and Mach number are denoted by Pe, Te, Ve and Me, respectively. The flow is locally subsonic at the convergent part of the nozzle sonic at the throat and supersonic at the exit of the nozzle.

1) Nozzle:

A nozzle is used to give designed to control the flow of the fluid (especially to increase the speed) at the exit point of the nozzle.

Nozzle is often a pipe with varying cross-sectional area, it can be used to direct or control the flow of fluid or modify the speed of the flow, they are used to modify the rate of flow in the pipe, direction, mass, shape, and pressure.

a basic diagram of 1D supersonic nozzle flow

Governing Equation:

There are 3 equations that govern the flow of the fluid in the nozzle continuity equation momentum equation and energy equation.

Non-Conservative Equation: In this type, a control volume is taken into consideration, where the control volume moves with the fluid and properties of the fluid are then taken into consideration.

CONTINUITY EQUATION:`-(delrho)/(delt) = - rho.(delv)/(delx)-rho.v.(dellnA)/(delx)-v.(delrho)/(delx)`

MOMENTUM EQUATION: `-(delv)/(delT) = -v.(delv)/(delx) - (1/rho).((delT)/(delx)+(t/rho).(delrho)/(delx))`

ENERGY EQUATION: `-(delT)/(delt) = -v.(delT)/(delx)- (1-rho).T.((delv)/(delx)+v.(dellnA)/(delx)`

CONSERVATIVE EQUATION: In this type, a control volume is taken into consideration where the control volume is fixed and the flow of the fluid in that control volume is studied.

CONTINUITY EQUATION: `- (del rhoA)/(delt) = -(delrhovA)/(delx)`

MOMENTUM EQUATION: `-(delrhoAv)/(delt) = - (del(rhoAv^2+(1/rho).rhoA))/(delx)+(1/rho).rho.(delA)/(delx)`

ENERGY EQUATION :`(del(rho.(T/(gamma-1)+gamma/2.v^2).A))/(delt)=-(del(rho.(T/(gamma-1)+gamma/2.v^2).v.A+rhoA.v))/(delx)`

MACCORMACK METHOD:

This method is used to solve the 1 dimensional quasi state supersonic nozzle flow . This method is a explicit finite difference technique which is second order accurate in both time and space.

In the case of computational fluid dynamics, the MacCormack method is used widely for discretization for numerical solutions for hyperbolic differential equations.

PREDICTOR STEP

In the predictor step, the terms in the right hand of the equation are changed to forward differencing method, the values of the left-hand side is calculated by time marching solution and for time-step \'t+dt\'.

This step is used will be used in the corrector step to rectify the errors and get correct results.

CORRECTOR STEP:

In this step, the terms in the right-hand side of the equation are changed to backward differencing method, the variables of the RHS are taken from the predictor step and are corrected in this step.

now, the predictor and corrected step are used to calculate the average time derivative which is used for the value of the variable for any given point for a particular time step.

clear all
close all
clc

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

%calculate 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 step
nt= 500;
dt=0.03;
tol=1e-3;
 
tic
%outer time loop
for k = 1:nt
  rho_old= rho;
  v_old = v;
  T_old = T;
  
  %predictor method
  
  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(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
    drho_dt_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)+drho_dt_p(j)*dt;
    v(j) = v(j)+dvdt_p(j)*dt;
    T(j) =T(j)+ dTdt_p(j)*dt;
  endfor
  %corrector method
    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(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
    drho_dt_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);
  endfor
%compute the average time derivative
drho_dt = 0.5*(drho_dt_p+drho_dt_c);
dvdt = 0.5*(dvdt_p+dvdt_c);
dTdt = 0.5*(dTdt_p+dTdt_c);

%final solution
for i = 2:n-1;
  rho(i) = rho_old(i)+drho_dt(i)*dt;
  v(i) = v_old(i)+dvdt(i)*dt;
  T(i) = T_old(i)+dTdt(i)*dt;
endfor
%Boundary Conditions
%inlet boundary conditons
v(1) = 2*v(2)-v(3);

%outlet boundary conditons
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);

%Mass flow rate
inflow = rho(1)*v(1)*a(1);
outflow = rho(n)*v(n)*a(n);
netmf = rho.*v.*a;

%Mach number
mach = max(v./sqrt(T));
% Error tolerance 
error_v(k) = max(v-v_old);
error_rho(k) = max(rho-rho_old);
error_T(k) = max(T-T_old);

if (error_v(k)< tol && error_rho(k) < tol
    &&error_T(k)<tol);
endif
endfor
time1 = toc;
%%plotting

%velocity plots
figure(1)
plot(x,v,\'b\');
xlabel(\'Distance in X-Direction\')
ylabel(\'Velocity\')
legend(\'Velocity in Non-Conservation\')
grid on
%density plot
figure(2)
plot(x,rho,\'r\')
xlabel(\'Distance in X-direction\')
ylabel(\'Density\')
legend(\'Density in Non-Conservation\')
grid on

%Temparature plot
figure(3)
plot(x,T,\'g\')
xlabel(\'Distance in x-direction\')
ylabel (\'Temprature\')
legend(\'Temprature in Non-Conservation\')
grid on

%Net mass flow rate plot
figure(4)
plot(x,netmf)
xlabel(\'Distance in x-direction\')
ylabel(\'net mass flow rate\')
legend(\'net mass flow rate in Non-Conservation\')
grid on

%Error in velocity with time
figure(5)
plot(error_v)
xlabel(\'Time\')
ylabel(\'Error in Velocity\')
legend(\'Error in Velocity\')
grid on

%Error in density with time
figure(6)
plot(error_rho)
xlabel(\'time step\')
ylabel(\'Error in Density\')
legend(\'Error in Density\')
grid on

%Error in Temparature with Time
figure(7)
plot(error_T)

xlabel(\'Time\')
ylabel(\'Error in Temparature\')
legend(\'Error in Temparature\')
grid on

Conservation code:

clear all
close all
clc

%input
n = 31;
x=linspace (0,3,n);
dx= x(2)-x(1);
%calculate initial profiles
rho = 1-0.3146*x;
T=1-0.2314*x;
gamma = 1.4;

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

%time step
nt = 500;

tol = 1e-6;

for i = 1:n
  deltat(i) = (0.5*(dx/((T(i)^0.5+v(i)))));
endfor
dt = 0.03;
%Derivatives 
U1= rho.*a;
U2=rho.*a.*v;
U3=rho.*(T/(gamma-1)+(gamma/2)*v.^2).*a;

tic
%outer time loop

for k=1:nt
  %storing old values
  U1_old = U1;
  U2_old = U2;
  U3_old = U3;
  
  rho_old = U1_old./a;
  v_old = U2_old./(rho.*a);
  T_old = (gamma-1)*((U3_old./U1_old)-(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 j = 2:n-1
    J2(j) = (1/gamma)*rho(j)*T(j)*(a(j+1)-a(j))/dx;
  endfor
  
  %predictor method
  for j = 2:n-1;
    
    %coninuity equation
    dU1dt_p(j) = -(F1(j+1)-F1(j))/dx;
    
    %momentum equation
    dU2dt_p(j) = -(F2(j+1)-F2(j))/dx +J2(j);
    
    %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;
  endfor
  
  %calculating variables
  rho = U1./a;
  v =U2./(rho.*a);
  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 j = 2:n-1
    J2(j) = (1/gamma)*rho(j)*T(j)*(a(j+1)-a(j))/dx;
  endfor
  
  %corrector method
  for j = 2:n-1
       %coninuity equation
    dU1dt_c(j) = -(F1(j)-F1(j-1))/dx;
    
    %momentum equation
    dU2dt_c(j) = -(F2(j)-F2(j-1))/dx + J2(j);
    
    %energy equation
    dU3dt_c(j) = -(F3(j)-F3(j-1))/dx;
  endfor
  %compute the 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 
  for i = 2:n-1;
    U1(i) = U1_old(i)+dU1dt(i)*dt;
    U2(i) = U2_old(i)+dU2dt(i)*dt;
    U3(i) = U3_old(i)+dU3dt(i)*dt;
  endfor
  
  %Boundary conditions
  
  %inlet
  U2(1) = 2*U2(2)-U2(3);
  
  %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);
  
  %recalculating variables
  rho = U1./a;
  v = U2./(rho.*a);
  T = (gamma-1)*((U3./U1)-(gamma/2)*v.^2);
  
  %Mass flow rate
  inflow = rho(1)*v(1)*a(1);
  outflow = rho(n)*v(n)*a(n);
  netmf = rho.*v.*a;
  
  %mach number 
  mach = v./sqrt(T);
  mach_number = max(mach);
  
  %error tolerance
  error_v(k) = max(v-v_old);
  error_rho(k) = max(rho-rho_old);
  error_T(k)= max(T-T_old);
  
  if (error_v(k)

PLOTS:

1) NON-CONSERVATIVE FORM:

fig 1- Velocity

fig 2- Density

fig 3- Tempertaure

fig 4- Net mass flow rate

mass flow near the throat is seen to be unstable as there are abrupt changes in the slope of the curve.

fig 5- Error in Velocity

fig 6- Error in Density

fig 7- Error in temperature

2) PLOTS FOR CONSERVATIVE FORM:

fig 1- Velocity

fig 2- Density

fig 3- Temperature

fig4 4- Net mass flow rate

The mass flow rate is more consistent in the throat area.

fig 5- Error in Velocity

fig 6- Error in Density

fig7 - Error in Temperature

RESULTS

GRID DEPENDENCE

after increasing the number of grid points it is observed that it has minimal effect on the result. When the grid point is doubled we reach close to the exact value but the margin is so small it is safe to say that the solution is independent of grid point.

The value of the time step is also required to be modified changing the grid point

n = 31 nt= 500

n = 61 nt= 1000

FASTEST METHOD

Time required for non-conservative form = 4.6734s

Time required for conservative form = 3.3565s

OBSERVATION

1) plots for velocity, temperature, density match with the exact plots

2) The net mass flow rate in the non-conservative form has more instabilities compared to conservative form, as in conservation flow the mass flow is a dependent variable whereas in the case of non-conservative form the mass flow is considered to be secondary and uses primitive variables to solve directly and closer to the exact solution

3) Solution varies in small margins after changing the grid points

4)Conservative form gives results faster than non-conservative form.

5)Mach number more than 1 suggests that the flow is supersonic.