Simulation of 1 Dimensional super-sonic nozzle flow
Objective
To write a code to solve the 1D supersonic nozzle flow equations using MacCormack method.
Problem Definition
The above image shows the different sections of the convergent divergent nozzle. The fluid flow is considered to be at steady state and Isentropic.
Governing Equations
The governing equations are solved in both conservative and Non-conservative form.
Non-conservative form
Mass conservation
`(delrho^')/(delt^')=-rho^'*(delv^')/(delx^')-rho^'*v^'*(dellnA')/(delx^')-v^'*(delrho^')/(delx^')`
Momentum conservation
`(delv^')/(delt^')=-v^'*(delv^')/(delx^')-1/gamma*((delT^')/(delx^')+T^'/rho^'*(delrho^')/(delx^'))`
Energy conservation
`(delT^')/(delt^')=-v^'*(delT^')/(delx^')-(gamma-1)*T^'((delv^')/(delx^')+v^'*(dellnA^')/(delx^'))`
Conservative form
Mass conservation
`(del(rho^'A^'))/(delt^')+(del(rho^'A^'v^'))/(delx^')=0`
Momentum conservation
`(del(rho^'A^'v^'))/(delt^')+(del[rho^'A^'v^('^2)+1/gamma*p^'A^'])/(delx^')=1/gamma*p^'*(delA^')/(delx^')`
Energy conservation
`(del[rho^'(e^'/(gamma-1)+gamma/2*v^('^2))*A^'])/(delt^')+(del[rho^'*(e^'/(gamma-1)+gamma/2*v^('^2))*v^'*A^'+p^'*A^'*v^'])/(delx^')=0`
solution vectors are defined as:
`u1=rho^'*A^', u2=rho^'*A^'*v^', u3=rho^'((T^')/(gamma-1)+gamma/2v^('^2))A^'`
Flux vectors are defines as:
`f1=u2`
`f2=(u2^2)/(u1)+1/gammau1(gamma-1)[(u3)/(u1)-gamma/2((u2)/(u1))^2]`
`f3=gamma(u2u3)/(u1)-(gamma(gamma-1))/2(u2^3)/(u1^2)`
The source term is defined by:
`j2=1/gammarho^'T^'(delA)/(delx)`
All the above equations are in dimensionless form.
MATLAB PROGRAM
Main program
% Quasi 1-D NOZZLE flow
clear all
close all
clc
% inputs
L=3; % length of the nozzle
n=31; % number of grid points
c=0.5; % Courant number
throat=(n-1)/2;
x=linspace(0,L,n); % Creating 1D grid
dx=x(2)-x(1); % grid space
gamma=1.4; % specific heat ratio
A=1+2.2*(x-1.5).^2; % Nozzle cross sectional area(function of x)
tic
nt=1400; % Number of time steps (Non-conservative form)
% Calling the function (Non-conservative form)
[rho,v,T,mach,mf,p]=non_conserv(x,dx,A,nt,gamma,c,n,throat);
sim_t_nc=toc;
tic
nt=550; % Number of time steps (Conservative form)
% Calling the function (Conservative form)
[rho2,v2,T2,mach_c,m_c,p_c]=conserv(A,x,nt,n,gamma,c,dx,throat);
sim_t_c=toc;
% Mass flow rate comparison between conservative and Non-conservative forms
figure(7)
plot(x,mf,'linewidth',2)
hold on
plot(x,m_c,'linewidth',2)
hold on
line([0 3],[0.579 0.579],'linewidth',2)
title('Mass flow rate comparison')
xlabel('Distance along Nozzle (Non-Dimensional)');
ylabel('Mass flow rate-(Non-Dimensional)');
legend('Non-conservative','conservative','exact solution')Non-Conservative form
1) Predictor step
function [d_rho_dt_p,d_v_dt_p,d_T_dt_p]=pred(rho,v,T,n,A,gamma,dx)
d_rho_dt_p=ones(1,31);
d_v_dt_p=ones(1,31);
d_T_dt_p=ones(1,31);
% Pre-allocating the array size
a=ones(1,31);
b=ones(1,31);
c=ones(1,31);
d=ones(1,31);
for i=2:n-1
a(i)=(v(i+1)-v(i))/dx;
b(i)=(log(A(i+1))-log(A(i)))/dx;
c(i)=(rho(i+1)-rho(i))/dx;
d(i)=(T(i+1)-T(i))/dx;
d_rho_dt_p(i)=-rho(i)*a(i)-rho(i)*v(i)*b(i)-v(i)*c(i);
d_v_dt_p(i)=-v(i)*a(i)-(1/gamma)*(d(i)+(T(i)/rho(i))*c(i));
d_T_dt_p(i)=-v(i)*d(i)-(gamma-1)*T(i)*(a(i)+v(i)*b(i));
end
2) Corrector step
function [drho_dt_c,dv_dt_c,dT_dt_c]=corr(rho,v,T,n,A,gamma,dx)
drho_dt_c=ones(1,31);
dv_dt_c=ones(1,31);
dT_dt_c=ones(1,31);
% Pre-allocating the array size
a=ones(1,31);
b=ones(1,31);
c=ones(1,31);
d=ones(1,31);
for i=2:n-1
a(i)=(v(i)-v(i-1))/dx;
b(i)=(log(A(i))-log(A(i-1)))/dx;
c(i)=(rho(i)-rho(i-1))/dx;
d(i)=(T(i)-T(i-1))/dx;
drho_dt_c(i)=-rho(i)*a(i)-rho(i)*v(i)*b(i)-v(i)*c(i);
dv_dt_c(i)=-v(i)*a(i)-(1/gamma)*(d(i)+(T(i)/rho(i))*c(i));
dT_dt_c(i)=-v(i)*d(i)-(gamma-1)*T(i)*(a(i)+v(i)*b(i));
end
3) Solving Non-Conservative equations
function[rho,v,T,mach,mf,p]=non_conserv(x,dx,A,nt,gamma,c,n,throat)
% Initial conditions
rho=1-(0.3146*x);
T=1-(0.2314*x);
v=(0.1+1.09*x).*T.^0.5;
% Pre-allocating the array size
p_throat=ones(1,nt);
rho_throat=ones(1,nt);
temp_throat=ones(1,nt);
mach_throat=ones(1,nt);
X=linspace(1,nt,nt);
% Iteration loop
for k=1:nt
% updating previous time_step values
rho_old=rho;
T_old=T;
v_old=v;
% time step calculation
sl=(sqrt(T)); % local velocity of sound
dt=c.*(dx./((sl)+v));
dt=min(dt);
% function for predictor step
[drho_dt_p,dv_dt_p,dT_dt_p]=pred(rho,v,T,n,A,gamma,dx);
for i=2:n-1
% next time step values
rho(i)=rho(i)+drho_dt_p(i)*dt;
v(i)=v(i)+dv_dt_p(i)*dt;
T(i)=T(i)+dT_dt_p(i)*dt;
end
% function for corrector step
[drho_dt_c,dv_dt_c,dT_dt_c]=corr(rho,v,T,n,A,gamma,dx);
% computing average values
d_rho_dt=0.5.*(drho_dt_p+drho_dt_c);
d_v_dt=0.5.*(dv_dt_p+dv_dt_c);
d_T_dt=0.5.*(dT_dt_p+dT_dt_c);
% Updating next time step values
for j=2:n-1
rho(j)=rho_old(j)+d_rho_dt(j)*dt;
v(j)=v_old(j)+d_v_dt(j)*dt;
T(j)=T_old(j)+d_T_dt(j)*dt;
end
% Boundary conditions (inlet)
v(1)=2*v(2)-v(3);
% Boundary conditions (exit)
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);
% Calculating state variables
mf=rho.*A.*v;
mach=v./(sqrt(T));
p=rho.*T;
% Calculating throat values
p_throat(k)=p(throat);
rho_throat(k)=rho(throat);
temp_throat(k)=T(throat);
mach_throat(k)=mach(throat);
figure(1)
if k==50
plot(x,mf,'color','r','linewidth',2)
hold on
elseif k==100
plot(x,mf,'color','g','linewidth',2)
hold on
elseif k==150
plot(x,mf,'color','b','linewidth',2)
hold on
elseif k==200
plot(x,mf,'color','bla','linewidth',2)
hold on
elseif k==700
plot(x,mf,'color','m','linewidth',2)
hold on
elseif k==1400
plot(x,mf,'y--o','linewidth',2)
legend('50-dt','100-dt','150-dt','200-dt','700-dt','1400-dt')
xlabel('Distance along Nozzle (Non-Dimensional)');
ylabel('Mass flow rate-(Non-Dimensional)')
title('Non-conservative form - Mass flow rate at different time steps')
hold on
end
end
figure(2)
subplot(2,1,1)
plot(X,rho_throat)
title('Non-conservative form - variation of density and temperature with time-steps at nozzle throat')
xlabel('No of time steps')
ylabel('Density')
subplot(2,1,2)
plot(X,temp_throat)
xlabel('No of time steps')
ylabel('Temperature')
figure(3)
subplot(2,1,1)
plot(X,p_throat)
title('Non-conservative form - variation of pressure and Mach number with time-steps at nozzle throat')
xlabel('No of time steps')
ylabel('Pressure')
subplot(2,1,2)
plot(X,mach_throat)
xlabel('No of time steps')
ylabel('Mach Number')
endConservative form
1) Forward Difference
function f1=fd_dx(f,dx,i)
% This function calculates the forward difference for a given function
f1=(f(i+1)-f(i))/dx;
end
2) Backward Difference
function f1=rd_dx(f,dx,i)
% This function calculates the rearward difference for a given function
f1=(f(i)-f(i-1))/dx;
end3) Solving Conservative equations
function [rho,v,T,mach_c,m_c,p_c]=conserv(A,x,nt,n,gamma,c,dx,throat)
% Pre-allocating the array size
rho=ones(1,n);
T=ones(1,n);
% Pre-allocating the array size
p_throat=ones(1,nt);
rho_throat=ones(1,nt);
temp_throat=ones(1,nt);
mach_throat=ones(1,nt);
X=linspace(1,nt,nt);
% Initial conditions
for i=1:n
if x(i)>=0&&x(i)<=0.5
rho(i)=1;
T(i)=1;
elseif x(i)>=0.5&&x(i)<=1.5
rho(i)=1-0.366*(x(i)-0.5);
T(i)=1-0.167*(x(i)-0.5);
elseif 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
% Defining velocity
v=0.59./(rho.*A);
% Defining Solution vector
u1=rho.*A;
u2=rho.*A.*v;
u3=rho.*A.*(T/(gamma-1)+gamma/2*v.^2);
% Iteration loop
for k=1:nt
% previous time step update
u1_old=u1;
u2_old=u2;
u3_old=u3;
% time step calculation
dt=min(c*dx./(sqrt(T)+v));
% Flux vectors
f1=u2;
f2=(u2.^2./u1)+((gamma-1)/gamma)*(u3-(gamma*u2.^2)./(2*u1));
f3=(gamma*u2.*u3./u1)-(gamma*(gamma-1)/2)*(u2.^3./u1.^2);
% Pre-allocating the array size
du1_dt_p=ones(1,n);
du2_dt_p=ones(1,n);
du3_dt_p=ones(1,n);
j2=ones(1,n);
% PREDICTOR STEP
for i=2:n-1
j2(i)=1/gamma*rho(i).*T(i)*fd_dx(A,dx,i);
du1_dt_p(i)=-fd_dx(f1,dx,i);
du2_dt_p(i)=-fd_dx(f2,dx,i)+j2(i);
du3_dt_p(i)=-fd_dx(f3,dx,i);
% solution update
u1(i)=u1(i)+du1_dt_p(i)*dt;
u2(i)=u2(i)+du2_dt_p(i)*dt;
u3(i)=u3(i)+du3_dt_p(i)*dt;
end
du1_dt_c=ones(1,n);
du2_dt_c=ones(1,n);
du3_dt_c=ones(1,n);
% updating primitive variables and flux vectors
rho=u1./A;
T=(gamma-1)*(u3./u1-(gamma/2)*(u2./u1).^2);
f1=u2;
f2=(u2.^2./u1)+((gamma-1)/gamma)*(u3-(gamma*u2.^2)./(2*u1));
f3=(gamma*u2.*u3./u1)-(gamma*(gamma-1)/2)*(u2.^3./u1.^2);
% CORRECTOR STEP
for j=2:n-1
j2(j)=1/gamma*rho(j).*T(j)*rd_dx(A,dx,j);
du1_dt_c(j)=-rd_dx(f1,dx,j);
du2_dt_c(j)=-rd_dx(f2,dx,j)+j2(j);
du3_dt_c(j)=-rd_dx(f3,dx,j);
end
% COMPUTING AVERAGE
du1_dt=0.5*(du1_dt_p+du1_dt_c);
du2_dt=0.5*(du2_dt_p+du2_dt_c);
du3_dt=0.5*(du3_dt_p+du3_dt_c);
% COMPUTING FINAL VALUES
for l=2:n-1
u1(l)=u1_old(l)+du1_dt(l)*dt;
u2(l)=u2_old(l)+du2_dt(l)*dt;
u3(l)=u3_old(l)+du3_dt(l)*dt;
end
% INLET BOUNDARY CONDITIONS
u2(1)=2*u2(2)-u2(3);
% OUTLET BOUNDARY CONDITION
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 of variables
rho=u1./A;
v=u2./u1;
T=(gamma-1)*(u3./u1-(gamma/2)*(u2./u1).^2);
% Calculating non-dimensional physical quantities
mach_c=v./sqrt(T);
m_c=rho.*A.*v;
p_c=rho.*T;
% Calculating throat values
p_throat(k)=p_c(throat);
rho_throat(k)=rho(throat);
temp_throat(k)=T(throat);
mach_throat(k)=mach_c(throat);
figure(4)
if k==50
plot(x,m_c,'color','r','linewidth',2)
hold on
elseif k==100
plot(x,m_c,'color','g','linewidth',2)
hold on
elseif k==150
plot(x,m_c,'color','b','linewidth',2)
hold on
elseif k==200
plot(x,m_c,'color','bla','linewidth',2)
hold on
elseif k==500
plot(x,m_c,'color','m','linewidth',2)
hold on
elseif k==550
plot(x,m_c,'y--o','linewidth',2)
legend('50-dt','100-dt','150-dt','200-dt','500-dt','600-dt')
xlabel('Distance along Nozzle (Non-Dimensional)');
ylabel('Mass flow rate-(Non-Dimensional)')
title('conservative form - Mass flow rate at different time steps')
hold on
end
end
figure(5)
subplot(2,1,1)
plot(X,rho_throat)
title('conservative form - variation of density and temperature with time-steps at nozzle throat')
xlabel('No of time steps')
ylabel('Density')
subplot(2,1,2)
plot(X,temp_throat)
xlabel('No of time steps')
ylabel('Temperature')
figure(6)
subplot(2,1,1)
plot(X,p_throat)
title('conservative form - variation of pressure and Mach number with time-steps at nozzle throat')
xlabel('No of time steps')
ylabel('Pressure')
subplot(2,1,2)
plot(X,mach_throat)
xlabel('No of time steps')
ylabel('Mach Number')
endResults
Minimum number of cycles for Convergence- (No of grid points=31)
The Non-conservative form takes more number of cycles for Convergence compared to Conservative form.
Mass Flow rate comparison between Conservative and Non-Conservative form (Grid points=31)
The Non-Conservative form has so many fluctuations compared to Conservative form in matching the exact solution.
Solution time
The time taken for calculation is the number of cycles to converge
31 Grid points
Solution time for conservative form is 0.572 sec
Solution time for Non-conservative form is 0.596 sec
61 Grid points
Solution time for conservative form is 1.100 sec
Solution time for Non-conservative form is 1.594 sec
The conservative form takes comparatively less number of cycles to converge but contains a lot of calculations while the Non-conservative form takes more number of cycles to converge but contains less calculations. However it depends on the code efficiency which depends on the structure of the code.
Grid dependence study
The fluctuations are less in the solution when the number of grid points is increased as seen in the above image.
Inference
The Conservative form of governing equations yield results closer to the exact solution compared to the Non-Conservative form of governing equations. A lot of fluctuations are also seen close to the inlet and outlet boundaries in the Non-conservative form of governing equation.