simulation of a quasi 1-d super-sonic sub-sonic nozzle flow using the macormack method

AIM:- To simulate a quasi 1-d flow through a converging - diverging nozzle by using the macormack method

quasi 1-d flow : A one dimensional as the name implies means that the flow is alog a given axis and the components of flow variables along other per[endicular axes are zero. that is in a 1-d flow along the x-axis the velocity component along the y and z axes will be zero. this type of flow rarely exists in real life and 2-d flows are most commonly encountered. but we can approximate the flow in a system in such a way that it resembles a 1-d flow. meaning we can assume the flow parameters along the other axes to be uniform along the respective axis. such type of flow is called a quasi one dimensional flow.

the convergent divergent nozle is connected to a reservoir where the pressure and temprature are `p_o`and `t_o`and are at stagnation points. the areaof the reservoir is assumed to be large making the velocity at inlet velocity zero.

flow at inlet : subsonic

flow at throat :sonic

flow at outlet :supersonic

here the area of the nozzle changes as a function of the distance qalong the nozzle,x, and therefore we can assume that the flow properties vary with x axis only and are uniform along the other axes thus makit it a quasi 1-d flow.

macormack method : it is an explicit finite difference technique which is second order accurate in both space and time. the plus point of macormack technique is that we can achieve second order accuracy without actually calculating the second order derivatives. this is done by calculating the derivative twice. first by forward differencing methos and then again by using the backward differencing method and finally taking their average to get the final value. the process can be completed in three steps

1. predictor step: the derivative of a quantity is calculated by using a forward differencing technique and solution is updated in following manner
   `(f_(i,j))^(t+deltat) = (f_(i,j))^(deltat) + ((deltaf)/(deltat))^t *nablat`
2. corrector step :the values obtained by the predictor steps are used to calculate the backward differencing derivative fot the quantity.
3. average solution: average is taken for the values obtained by the predictor and corrector steps and we get the corrected value of the quantity at time `t+nablat`((delf)/(delt))_(avg) = 1/2*[(((delf)/(delt))_(i,j))^t +(((delf)/(delt))_(i,j))^(t+nablat) ]`final solution is updated by `(f_(i,j))^(t+deltat) = (f_(i,j))^(deltat) + (((deltaf)/(deltat))_(avg))^t *nablat`

because a two step predictor -corrector sequence is used with forward differencing on predictor and backward differencing on the corrector we get a second order accurate solution.

in macormack method the differencing techniques used predictor and corrector steps are not a rule of thumb. same order accuracy will be obtained even if backward differencing is done on the predictor step and forward differencing is done in corrector step.

without the physically proper implementation of boundary conditions and their numerically proper representation we cannot obtain a proper numerical solution.point 1 is essentially in the reservoir and it represents an inflow boundary, with flow coming from the reservoir and into the nozzle and point n is an outflow boundary, with flow leaving the nozzle at the nozzle exit. the method of characterstics tells us that two conditions must be specified at the boundary and other dependent flow variable must be allowed to float at the boundary

time step calculations: since the macormack method is an explicit approach, its stability depends greatly on the time step used for the calculations. we can use a CFL criterion for calculating a stable time step.

`nablat = (nablax)/(a+v)`

by using this criteria we can determine the optimal time step for each local grid point. but this creates a problem, the resulting flow field will be in a type of artificial \"time warp\", and we would get non physical data and the solution will bmnot be reliable. so to solve this problem we calculate the time step for all the points and choose the minimum value from the obtained values.

Non Conservative approach : the control volume for which flow is not stationary.

the system can be represented by following set of non linear partial differential equations:

1. continuity equation :`(delrho)/(delt) = -rho*(delv)/(delx) - rhovdel(lna)/(delx) - v(delrho)/(delx)`
2. momentum equation: `-v*(delv)/(delx) - (1/gamma)*[(delT)/(delx) + (T/rho)(delp)/(delx))]`
   
3. energy equation: `(delT)/(delt) = -v(delT)/(delx) - (gamma-1)*T[(delv)/(delx) +v(del(lna))/(delx)]`

the above equations are not represented in their standard forms instead they are written in non dimensional form. it means that the variables in the equation does not represents their respective quantities but they represents the ratios of said quantities with their values at stagnation point. the non dimensional values vary between 0 and 1 which is an aesthetic advantage while presenting the result. the simulation can be performed using the standard equations also without any change in the results.

initial conditions:

• area `A = 1+2.2(x/l -1.5)^2`
• density `rho = 1-0.3146x`
• temprature `T = 1-0.2314x`
• velocity = `v= (0.1+1.09x)T^(1/2)`

Conservative approach: the control volume is assumed to be stationary

the system is represented by following set of non linear partial differential equations:

1. continuity equation = `(del(rhoav))/(delt) = (del(rhoav))/(delx)`
2. momentum equation = `(del(rhoav))/(delt) + (del/(delx))[rhoav^2 +(1/gamma)pa] = (1/gamma)p((dela)/(delx))`
3. energy equation`del/(delt)[rhoa(e/(gamma-1) + gamma/2v^2) ] + del/(delx)[rhoav(e/(gamma-1)+gamma/2v^2+pav )] = 0`

the above stated equations when simulated will give garbage values after about 300 steps. this is termed as a \'lack of purity\' so to avoid this we can formulate the equations in another manner. by assuming some new variables which are functions of the primitive variables we can simplify the equations and stablise the solution. these terms are called the solution vectors,flux vectors, and source terms. these flux terms can be solved and the values of the primitives can be calculated from the obtained flux values.

variables

u = solution vectors

f = flux vectors

j = source terms

• `u_1 = rhoa`
• `u_2 = rhoav`
• `u_3 = rhoa(e/gamma-1 + gamma/2v^2)`
• `f_1 = u_1`
• `f_2 = u_2^2/u_1 + (gamma-1)/gamma(u_3 -gamma/2u_2^2/u_1)`
• `f_3 = gamma(u_3u_2)/u_1 - gamma((gamma-1))/2(u_2^3/u_1^2)`
• `j = (gamma-1)/gamma(u_3-gamma/2u_2^2/u_1)del/(delx)(lna)`

by replacing the primitives with these variables we can make the equation numerically \'pure\'

1. continuityequation:`(delu_1)/(delt) = -(delf_1)/(delx)`
2. momentum equation:`(delu_2)/(delt) = -(delf_2)/(delx) + j`
3. energy equation:`(delu_3)/(delt) = -(delf_3)/(delx)`

initial conditions:

for `0<=x<=0.5`

• `rho = 1.0`
• `T = 1.0`

for `0.5<=x<=1.5`

• `rho = 1.0-0.366(x-0.5)`
• `T = 1.0-0.167(x-0.5)`

for`1.5<=x<=3.5`

• `rho = 0.634-0.3879(x-1.5)`
• `T = 0.833-0.3507(x-1.5)`

assuming a constant mass flow rate that is `u_2 = rhoav` close to the analytic value of 0.59. therefor `v = 0.59/(rhoa)`

the conservative form equations are theoretically identical to those which are derived in terms of soltion and flux vectors. so logically these both equations should give same results but the equation in terms of primitives blows and the one with generic terms the solution is stable!! mathematically this should not happen. the reason which i believe this happens is because the truncation error grows exponentially in the primitive equations.

Boundary conditions : the boundary conditions does not depends on the approach of the solution that is wether we are employing non conservative or conservative approach. basically two boundaries are fixed at the subsonic inflow boundary and one condition is allowed to float and at the supersonic boundary all properties are allowed to float

• `v_n = 2v_(n-1) -v_(n-2)`
• `rho_n = 2rho_(n-1) -rho_(n-2)`
• `T_n = 2T_(n-1) -T_(n-2)`

it can be observed from the initial conditions that in case of non conservative approach the boundary is defined quite loosely than compared to the conservative approach, so it should be obvious that the conservative approach will give better solution with less errors and closer to the analytical solution.

code:

main code:

clear all
close all
clc

% inputs
n = 31;
x = linspace(0,3,n);
dx = x(2) - x(1);
a = 1+2.2*(x-1.5).^2;
th = (n+1)/2;
gamma = 1.4;

c = 0.5;
nt = 1400;

%function for solving the problem by non conservative method
tic
[v_n,rho_n,t_n,m_n,mach_n,p_n] = non_conservative_fun(a,gamma,n,dx,nt,x,c);
t1 = toc;

%function for solving the problem by conservative method
tic
[v_c,rho_c,t_c,u2,mach_c,p_c] = conservative_fun(a,gamma,n,dx,nt,x,c);
t2 = toc;

%plotting 
figure(3)
%comparison of mass flow variation obtained with the conservative and non conservative forms
plot(x,m_n)
hold on
plot(x,u2)
title(\'conservative vs non conservative mass flow rates\')
legend(\'non-conservative\',\'conservative\')
xlabel(\'length\')
ylabel(\'mass flow rate\')

%plots the values obtained by non conservative equations
figure(4)
title(\'by non conservative approach\')
subplot(4,1,1)
plot(x,mach_n)
title(\'variation of mach number\')
xlabel(\'length\')
ylabel(\'mach number\')

subplot(4,1,2)
plot(x,rho_n)
title(\'variation of density\')
xlabel(\'length\')
ylabel(\'density\')

subplot(4,1,3)
plot(x,p_n)
title(\'variation of pressure\')
xlabel(\'length\')
ylabel(\'pressure\')

subplot(4,1,4)
plot(x,t_n)
title(\'variation of temprature\')
xlabel(\'length\')
ylabel(\'temprature\')

%plots the values obtained by solving the conservative equations
figure(5)
title(\'by conservative approach\')
subplot(4,1,1)
plot(x,mach_c)
title(\'variation of mach number\')
xlabel(\'length\')
ylabel(\'mach number\')

subplot(4,1,2)
plot(x,rho_c)
title(\'variation of density\')
xlabel(\'length\')
ylabel(\'density\')

subplot(4,1,3)
plot(x,p_c)
title(\'variation of pressure\')
xlabel(\'length\')
ylabel(\'pressure\')

subplot(4,1,4)
plot(x,t_c)
title(\'variation of temprature\')
xlabel(\'length\')
ylabel(\'temprature\')

the code defines the domain and then calls the respective functions to solve by conservative and non conservative approaches. the values thus obtained are then plotted in form of a table to better understand the trend which is being followed by the properties. another plot is obtained which shows the comparison between the conservative and non conservative approaches.the time elapsed in running of each function is stored for comparison.

function for non conservative approach:

function [v_n,rho_n,t_n,m_n,mach_n,p_n]  = non_conservative_fun(a,gamma,n,dx,nt,x,c)
%initial conditions

  rho_n = 1-0.3146*x;
  t_n = 1-0.2314*x;
  v_n = (0.1 + 1.09.*x).*t_n.^(0.5);
  

  
  %outer time loop
  for k = 1:nt
    
    rho_old = rho_n;
    t_old = t_n;
    v_old = v_n;

  
    %defining time step using the CFL criteria
  %the time step at a given grid point i is given by
  %dt = (c*dx)/(a(i) + v(i))
  A = sqrt(t_n);
  dt = c*dx./(A+v_n);
  dt = min(dt);

  
  %predictor step
  %using forward differencing to calculate derivatives
  for j = 2:n-1
  %defining the temporary variables
  dvdx(j) = (v_n(j+1) - v_n(j))/dx;          %derivative of velocity w.r.t length
  drhodx(j) = (rho_n(j+1) - rho_n(j))/dx;    %derivative of density w.r.t length
  dtdx(j) = (t_n(j+1) - t_n(j))/dx;            %derivative of temprature w.r.t length
  dlogadx(j) = (log(a(j+1)) - log(a(j)))/dx; %derivative of log of area w.r.t length
  
  %continuity equation
  drhodt_p(j) = -rho_n(j)*dvdx(j) - rho_n(j)*v_n(j)*dlogadx(j) - v_n(j)*drhodx(j);
  %momentum equation
  dvdt_p(j) = -v_n(j)*dvdx(j) - (1/gamma)*(dtdx(j) + (t_n(j)/rho_n(j))*drhodx(j));
  %energy equation
  dtdt_p(j) = -v_n(j)*dtdx(j) - (gamma-1)*t_n(j)*(dvdx(j) + v_n(j)*dlogadx(j));
  
  %solution update
  rho_n(j) = rho_n(j) + drhodt_p(j)*dt;
  v_n(j) = v_n(j) + dvdt_p(j)*dt;
  t_n(j) = t_n(j) + dtdt_p(j)*dt;
  
end

%corrector step
  %using backward differencing to calculate derivatives
  for j = 2:n-1
  %defining the temporary variables
  dvdx(j) = (v_n(j) - v_n(j-1))/dx;          %derivative of velocity w.r.t length
  drhodx(j) = (rho_n(j) - rho_n(j-1))/dx;    %derivative of density w.r.t length
  dtdx(j) = (t_n(j) - t_n(j-1))/dx;            %derivative of temprature w.r.t length
  dlogadx(j) = (log(a(j)) - log(a(j-1)))/dx; %derivative of log of area w.r.t length
  
  %continuity equation
  drhodt_c(j) = -rho_n(j)*dvdx(j) - rho_n(j)*v_n(j)*dlogadx(j) - v_n(j)*drhodx(j);
  %momentum equation
  dvdt_c(j) = -v_n(j)*dvdx(j) - (1/gamma)*(dtdx(j) + (t_n(j)/rho_n(j))*drhodx(j));
  %energy equation
  dtdt_c(j) = -v_n(j)*dtdx(j) - (gamma-1)*t_n(j)*(dvdx(j) + v_n(j)*dlogadx(j));
end

%computing average values
  drhodt = 0.5*(drhodt_p + drhodt_c);
  dvdt = 0.5*(dvdt_p + dvdt_c);
  dtdt = 0.5*(dtdt_p + dtdt_c);
  
%final solution update
  for i = 2:n-1
  rho_n(i) = rho_old(i) + drhodt(i)*dt;
  v_n(i) = v_old(i) + dvdt(i)*dt;
  t_n(i) = t_old(i) + dtdt(i)*dt;
  end  

  %boundary conditions
  v_n(1) = 2*v_n(2) - v_n(3);
  %outlet conditions
  rho_n(n) = 2*rho_n(n-1) - rho_n(n-2);
  v_n(n) = 2*v_n(n-1) - v_n(n-2);
  t_n(n) = 2*t_n(n-1) - t_n(n-2);
  
  %calculations
  %non dimensionless mass flow rate
  m_n = rho_n.*a.*v_n;
  %mach number
  mach_n = v_n./(sqrt(t_n));
  %non dimensionless pressure
  p_n = rho_n.*t_n;
  
  
  %plot(x,v_n)
  %pause(0.03)
  %plotting
  figure(1)
  grid on
  hold on
  if k==50
    plot(x,m_n,\'r\')
    else if k==100
    plot(x,m_n,\'b\')
    else if k==150
    plot(x,m_n,\'g\')
   else if k==200
    plot(x,m_n,\'y\')
   else if k==700
    plot(x,m_n,\'k\')
   else if k==1400
    plot(x,m_n,\'c\')
   title(\'distribution of dimensionless mass flow rate along the length of the nozzle calculated in non conservative form\')
    legend(\'50dt\',\'100dt\',\'150\',\'200dt\',\'700dt\',\'1400dt\')
    xlabel(\'x/l\')
    ylabel(\'dimensionless mass flow rate\')
   end
  end
 end
end
end
end
end
endfunction

the function takes in the necessary inputs from the main code and defines the initial conditions. the outer time loop which corresponds to the time marching of the solution executes the entire code till specified number of times.time step for the solution is calculated inside the time loop so that with each iteration the time step gets updated with the other values. the code consists of standard steps for maccormack method, being the predictor, corrector,solution update and the average steps. to simplify the equations which are being solved some temporary variables are created so that it is easy to go through the code for debugging purposes. boundary conditions are specified to finish the problem statement. after the primitive values are calculated relevent quantities such as mass flow rate , mach number and pressure are calculated. a graph is then plotted against the length of the nozzle and the mass flow rate at different number of time the time marching is performed.

code for conservative approach:

function [v_c,rho_c,t_c,u2,mach_c,p_c] = conservative_fun(a,gamma,n,dx,nt,x,c)

%initial conditions
for i = 1:n
  if(0<=x(i)&&x(i)<=0.5)
    rho_c(i) = 1;
    t_c(i) = 1;
  end
  
  if (0.5<=x(i)&&x(i)<=1.5)
    rho_c(i) = 1-0.366*(x(i)-0.5);
    t_c(i) = 1-0.167*(x(i)-0.5);
  end
  
  if (1.5<=x(i)&&x(i)<=3.5)
    rho_c(i) = 0.634-0.3879*(x(i)-1.5);
    t_c(i) = 0.833-0.3507*(x(i)-1.5);
  end
  
endfor

%the conservative form of the equation is not solved for the primitive variables
%instead primitives are replaced by variables which are function of them
%these variables are solved for and then the primitive values are calculated from them
 
v_c = 0.59./(rho_c.*a);
u1 = rho_c.*a;
u2 = u1.*v_c;
u3 = u1.*((t_c/(gamma-1))+ (gamma/2).*(v_c.^2));

  %outer time loop
  for k = 1:1400
  %defining time step using the CFL criteria
  %the time step at a given grid point i is given by
  %dt = (c*dx)/(a(i) + v(i))
  A = sqrt(t_c);
  dt = c*dx./(A+v_c);
  dt = min(dt);
  
  u1_old = u1;
  u2_old = u2;
  u3_old = u3;
  
  %determining the flux values
  f1 = u2;
  term1 = (u2.^2)./u1;
  term2 = u3-((gamma/2)*term1);
  f2 = term1+((gamma-1)/gamma)*term2;
  term3 = (u2.*u3)./u1;
  term4 = (u2.^3)./(u1.^2);
  f3 = gamma.*term3-((gamma*(gamma-1))/2).*term4;
  
  %predictor method
  for j = 2:n-1
    du1dt_p(j) = -(f1(j+1)-f1(j))/dx;
    du2dt_p(j) = (-(f2(j+1)-f2(j)))/dx +(1/gamma)*rho_c(j)*t_c(j)*((a(j+1)-a(j))/dx);
    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
  
  %determining the predicted primitive values
  rho_c = u1./a;
  v_c = u2./u1;
  t_c = (gamma-1)*((u3./u1) - (gamma/2)*(v_c.^2));
  %determining the flux values
  f1 = u2;
  term1 = (u2.^2)./u1;
  term2 = u3-((gamma/2)*term1);
  f2 = term1+((gamma-1)/gamma)*term2;
  term3 = (u2.*u3)./u1;
  term4 = (u2.^3)./(u1.^2);
  f3 = gamma.*term3-((gamma*(gamma-1))/2).*term4;
  
  %corrector step
  for j = 2:n-1
    du1dt_c(j) = -(f1(j)-f1(j-1))/dx;
    du2dt_c(j) = (-(f2(j)-f2(j-1)))/dx +(1/gamma)*rho_c(j)*t_c(j)*((a(j)-a(j-1))/dx);
    du3dt_c(j) = -(f3(j)-f3(j-1))/dx;  
  end
  
  %computing average
  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 l = 2:n-1
    u1(l) = u1_old(l) + du1dt(l)*dt;
    u2(l) = u2_old(l) + du2dt(l)*dt;
    u3(l) = u3_old(l) + du3dt(l)*dt;
  end

  %inlet boundary conditions
  u1(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));
  
  %outlet boundary conditions
  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);
  
  %determining the primitive variables
  rho_c = u1./a;
  v_c = u2./u1;
  t_c = (gamma-1)*((u3./u1) - (gamma/2)*(v_c.^2));
  mach_c = v_c./(t_c.^(0.5));
  p_c = rho_c.*t_c;
  
  %plotting
  %plots the variation of mass flow rate at different time steps 
  figure(2)
  grid on
  hold on
    if k==50
    plot(x,u2,\'r\')
    else if k==100
    plot(x,u2,\'b\')
    else if k==150
    plot(x,u2,\'g\')
   else if k==200
    plot(x,u2,\'y\')
   else if k==700
    plot(x,u2,\'k\')
   else if k==1400
    plot(x,u2,\'c\')
   title(\'distribution of dimensionless mass flow rate along the length of the nozzle calculated in conservative form\')
    legend(\'50dt\',\'100dt\',\'150\',\'200dt\',\'700dt\',\'1400dt\')
    xlabel(\'x/l\')
    ylabel(\'dimensionless mass flow rate\')
   end
  end
 end
end
end
end
end
endfunction

the code takes in necessary inputs from the main code and initial conditions are calculated. as discussed above it is evident that the initial conditions are more defined. the primitive equations are then converted into generic equations in terms of solution vector, flux vectors and source terms. these generic equations are then fed into he macormack method to get the values for solution vectors. final values of primitive variables are then calculated from these solution vectors. a plot is generated between mass flow rate at different time steps and the dimensionless length of the nozzle.

plots:

1. plot between mass flow rate and time steps obtained by non conservative approach

2.plot between mass flow rate and time steps obtained by conservative approach

3.comparison between conservative and non conservative solutions:

4.data obtained from non conservative approach:

5.data obtained by conservative approach:

grid dependency:

the values for flow parameters like mach number, flow rate,density and pressure is calculated for both conservative and non conservative at the throat of the nozzle and compared at different grid points.

for non conservative solution:

	grid dependency test
grid no.	flow rate	mach no.	density	pressure
31	0.58377	0.99937	0.63869	0.53424
61	0.58255	0.9997	0.63756	0.53261

for conservative solution:

	grid dependency test
grid no.	flow rate	mach no.	density	pressure
31	0.58522	0.98173	0.65039	0.54635
61	0.58287	0.99935	0.63821	0.53302

the difference between the respective flow parameters at the throat of the nozzle is not very large, so it is safe to assume that the mrthod is independent of the number of the grid points for which the solution is obtained

conclusions:

• macormack method is a highly accurate method.
• the solution becomes stable about 500 to 700 iterations.
• the solution obtained by the conservative approach is more accurate than the non conservative one. partly due to initial conditions being more efined in this case.
• the trend in the variation of mach number, density, pressure and temprature with the length of the nozzle is as expected.
• conservative approach is more accurate but it also requires more computational effort
• the solution is grid independent