Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
Problem Statement: Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method.
- To Figure out the minimum number of cycles for which the simulation should be run in order for convergence.
- To Compare the normalized mass flow rate between the conservative forms and non-conservative forms and make suitable conclusions.
- To find the faster method.
1. Topic: Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method.
2. Theory: The present problem deals with quasi-one-dimensional flow through a convergent-divergent nozzle.
a) Physical Description: The flow is considered to be steady, isentropic flow through a convergent-divergent nozzle as shown in figure below. The flow at inlet comes fram a reservoir where pressure, density and Temperature can be taken as P0, ρ0 and T0 respectively.
The flow is locally subsonic in the convergent section of the nozzle, sonic at throat (minimum Area) and supersonic at divergent section. An assumption has been made that the flow properties vary only with x.
Governing Equations: These are basically the mathematical statements of three fundamental physics principles :
1. THE CONTINUITY EQUATION: It is based on the conservation of mass principle. It can be expressed as:
ρ1V1A1=ρ2V2A2
2. THE MOMENTUM EQUATION: It is based on Newton\'s Second law and can be expressed as:
`p1*A1 + ρ1*V1^2*A1 + ∫_(A1 to A2) pdA = p2*A2 + ρ2*V2^2*A2`
3. THE ENERGY EQUATION: It is based on the principle the energy is conserved and can be expressed as:
h1+(V1^2)/2 = h2+(V2^2)/2
These governing equations can be derived into two forms i.e. conservative and Non-conservative forms.
The technique used to perform CFD solution of Subsonic- Supersonic Isentropic Nozzle is MACCORMACK\'S Technique.
Maccormack\'s Technique: It is an explicit finite-difference technique that is second-order accurate in both space and time. It is a two step technique with predictor and corrector step.
3. Program:
a) Non-conservation: It is derived for the infinitesimally small fluid element is moving with the flow. In this method, the governing non-conservative dimensionless equations derived are:
Continuity Equation:
`∂ρ/∂t= -ρ*∂V/∂x - ρ*V*(∂(ln〖A〗 )/∂x – V*∂ρ/∂x )`
Momentum Equation:
`∂V/∂t= -V*∂V/∂x – 1/γ( ∂T /∂x + T/ρ*∂ρ/∂x )`
Energy Equation:
`∂T/∂t= -V*∂T/∂x – ( γ-1 )*T*( ∂V/∂x + V*∂(ln〖A〗)/∂x )`
Calculation of Time step: It is done by the Courant-Friedrichs-Lowry criterion for 1-d flow, which is given by;
`∆t = C*∆x/(a+V)`
It is calculated at all grid points and the minimum value is finally chosen.
Boundary Conditions:
Subsonic Inflow boundary:
ρ1=1
T1=1
V1=2*V2-V3
Supersonic outflow Boundary:
Vn=2*V(n-1)-V(n-2)
ρn= 2*ρ(n-1)- ρ(n-2)
Tn= 2*T(n-1)- T(n-2)
Nozzle Shape and Initial conditions:
Nozzle Area is specified and held fixed, independent of time.
A=1+2.2*(x-1.5)^2
Also as the flow expands through the nozzle; density and Temp. decreases and V increases and the initial conditions are given by;
ρ=1-0.3146*x
T=1-0.2314*x
V=(0.1+1.09*x)*T^(1/2)
at t=0;
function [rho,v,T,m,p,rhot1,vt1,Tt1]= non_conver(n,x,dx,A,g,nt,trt,C)
%calculate initial profiles
rho=1-0.3146*x;
T=1-0.2314*x;
v=(0.1+1.09*x).*T.^0.5;
dt=C*min(dx./((T.^0.5)+v));
%Outer Time Loop
for k=1:nt
rho_old=rho;
T_old=T;
v_old=v;
A_old=A;
%Predictor Time step
for j=2:n-1
dvdx=(v(j+1)-v(j))/dx;
drhodx=(rho(j+1)-rho(j))/dx;
dadx=(log(A(j+1))-log(A(j)))/dx;
dTdx=(T(j+1)-T(j))/dx;
%Governing Equations
drhodt(j)=(-rho(j)*dvdx)-(rho(j)*v(j)*dadx)-v(j)*drhodx;
dvdt(j)=(-v(j)*dvdx)-(dTdx+T(j)*drhodx/rho(j))/g;
dTdt(j)=(-v(j)*dTdx)-((g-1)*T(j)*(dvdx+v(j)*dadx));
%Updation
rho(j)=rho_old(j)+drhodt(j)*dt;
v(j)=v_old(j)+dvdt(j)*dt;
T(j)=T_old(j)+dTdt(j)*dt;
end
for(j=2:n-1)
%Corrector
dvdx=(v(j)-v(j-1))/dx;
drhodx=(rho(j)-rho(j-1))/dx;
dadx=(log(A(j))-log(A(j-1)))/dx;
dTdx=(T(j)-T(j-1))/dx;
%Governing Equations
drhodt_c(j)=(-rho(j)*dvdx)-(rho(j)*v(j)*dadx)-v(j)*drhodx;
dvdt_c(j)=(-v(j)*dvdx)-(dTdx+T(j)*drhodx/rho(j))/g;
dTdt_c(j)=(-v(j)*dTdx)-((g-1)*T(j)*(dvdx+v(j)*dadx));
end
drhodt_a=(drhodt+drhodt_c)./2;
dvdt_a=(dvdt+dvdt_c)./2;
dTdt_a=(dTdt+dTdt_c)./2;
for j=2:n-1
rho(j)=rho_old(j)+drhodt_a(j)*dt;
v(j)=v_old(j)+dvdt_a(j)*dt;
T(j)=T_old(j)+dTdt_a(j)*dt;
end
%Boundary Conditions
v(1)=2*v(2)-v(3);
rho(1)=1;
T(1)=1;
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
m= rho.*A.*v;
%Pressure
p=rho.*T;
%Throat Values
rhot1(k)=rho(trt);
vt1(k)=v(trt);
Tt1(k)=T(trt);
%Graph
figure(1)
if k==100
plot(x,m,\'r\',\'linewidth\',1.5);
hold on
elseif k==200
plot(x,m,\'b\',\'linewidth\',1.5);
hold on
elseif k==300
plot(x,m,\'g\',\'linewidth\',1.5);
hold on
elseif k==400
plot(x,m,\'y\',\'linewidth\',1.5);
hold on
elseif k==500
plot(x,m,\'k\',\'linewidth\',1.5);
hold on
elseif k==600
plot(x,m,\'c\',\'linewidth\',1.5);
hold on
elseif k==700
plot(x,m,\'m\',\'linewidth\',1.5);
hold on
elseif k==800
plot(x,m,\'*r\',\'linewidth\',1.5);
hold on
elseif k==900
plot(x,m,\'*b\',\'linewidth\',1.5);
hold on
elseif k==1000
plot(x,m,\'*g\',\'linewidth\',1.5);
grid on
xlabel(\'X\');
ylabel(\'Mass flow rate\');
title(\'Non-conversative form mass flow rate along with x \',\'color\',\'k\');
legend({\'100 time steps\',\'200 time step\',\'300 time step\',\'400 time step\',\'500 time step\',\'600 time step\',\'800 time step\',\'900 time step\',\'1000 time step\'});
end
end
endb) Conservative Form: It is derived for the infinitesimally small fluid element that is fixed in space. In this method, the governing non-conservative dimensionless equations derived are:
Continuity Equation:
`(∂(ρ*A))/∂t+ (∂(ρ*A*v))/∂x = 0`
Momentum Equation:
`(∂(ρ*A*v))/∂t+ (∂[ρ*A*v^2+ 1/γ*ρ*A])/∂x = 1/γ*ρ*(∂A )/∂x`
Energy Equation:
`(∂(ρ*(e/(γ-1)+ γ/2*v^2 )*A))/∂t + (∂(ρ*(e/(γ-1)+ γ/2*v^2 )*v*A+ρ*A*v))/∂x = 0`
As these equations are not in generic form; they can be expressed as such using solution vectors U, flux vectors F, and source term J.
U1= ρ*A
U2= ρ*A*v
U3= ρ*(e/(γ-1)+ γ/2*v^2 )*A
F1= ρ*A*v
F2= ρ*A*v^2+ 1/γ*ρ*A
F3= ρ*(e/(γ-1)+ γ/2*v^2 )*v*A+ρ*A*v
J2= 1/γ*ρ*(∂A)/∂x
So the equations become;
`∂U1/∂t= -∂F1/∂x`
`∂U2/∂t= -∂F2/∂x +J2`
`∂U3/∂t= -∂F3/∂x`
Calculation of Time step: It is done by the Courant-Friedrichs-Lowry criterion for 1-d flow, which is given by;
`∆t = C*∆x/(a+V)`
It is calculated at all grid points and the minimum value is finally chosen.
Boundary Conditions:
Subsonic Inflow boundary:
U1(i=1)= ρ*A (at i=1)
U2(i=1)=2*U2(i=2)-U2(i=3)
U3(i=1)=U1(T/(γ-1)+ γ/2*v^2) (at i=1)
Supersonic outflow Boundary:
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)
Nozzle Shape and Initial conditions:
Nozzle Area is specified and held fixed, independent of time.
A=1+2.2*(x-1.5)^2
Also as the flow expands through the nozzle; density and Temp. decreases and V increases and the initial conditions are given by;
For x=0 to 0.5
ρ=1
T=1
For x=0.5 to 1.5
ρ=1-0.366*(x-0.5)
T=1-0.167*(x-0.5)
For x=1.5 to 3.5
ρ=0.634-0.3879*(x-1.5)
T=0.833-0.3507*(x-1.5)
V=U2(i=1)/(ρ*A)
at t=0;
function [rho,v,T,m,p,rhot2,vt2,Tt2]= conserve(n,x,dx,A,g,nt,trt,C)
%calculate initial profiles
for i=1:n
if (x(i)<=0.5)
rho(i)=1;
T(i)=1;
elseif (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)<=3.5)
rho(i)=0.634-0.3879*(x(i)-1.5);
T(i)=0.833-0.3507*(x(i)-1.5);
end
end
v=0.59./(rho.*A);
%CFL
dt=C*min(dx./((T.^0.5)+v));
%Solutions vectors
u1=rho.*A;
u2=rho.*A.*v;
u3=u1.*(T./(g-1)+((g/2).*(v.^2)));
%Outer Time Loop
for j=1:nt
u1_old=u1;
u2_old=u2;
u3_old=u3;
%Flux Terms
f1=u2;
j2=((g-1)/g)*(u3-(g/2*((u2.^2)./u1)));
f2=(((u2.^2)./u1)+j2);
f3=(((g*u2.*u3)./u1)-((g*(g-1)*(u2.^3))./(2*u1.^2)));
%Inner Loop
for k=2:n-1
%Predictor Step
dladx(k)=(log(A(k+1))-log(A(k)))/dx;
dadx(k)=(A(k+1)-A(k))/dx;
j_term(k)= (1/g)*(rho(k)*T(k))*dadx(k);
%Continuity
du1dt_p(k)=-(f1(k+1)-f1(k))/dx;;
%Momentum
du2dt_p(k)=-(f2(k+1)-f2(k))/dx+ j_term(k);
%Energy
du3dt_p(k)=-(f3(k+1)-f3(k))/dx;
u1(k)=u1(k)+du1dt_p(k)*dt;
u2(k)=u2(k)+du2dt_p(k)*dt;
u3(k)=u3(k)+du3dt_p(k)*dt;
end
rho=u1./A;
v=(u2./u1);
T=(g-1)*((u3./u1)-((g/2)*((u2./u1).^2)));
f1=u2;
j2=((g-1)/g)*(u3-(g/2*(u2.^2)./u1));
f2=(((u2.^2)./u1)+j2);
f3=(((g*u2.*u3)./u1)-((g*(g-1)*(u2.^3))./(2*u1.^2)));
%Forward Diffreencing
for k=2:n-1
%Predictor step
df1dx(k)= (f1(k)-f1(k-1))/dx;
df2dx(k)= (f2(k)-f2(k-1))/dx;
df3dx(k)= (f3(k)-f3(k-1))/dx;
dladx(k)=(log(A(k))-log(A(k-1)))/dx;
dadx(k)=(A(k)-A(k-1))/dx;
j_term2(k)= (1/g)*(rho(k)*T(k))*dadx(k);
%Continuity
du1dt_c(k)=-df1dx(k);
%Momentum
du2dt_c(k)=-df2dx(k)+ j_term2(k);
%Energy
du3dt_c(k)=-df3dx(k);
end
du1dt_av=(du1dt_p+du1dt_c)*0.50;
du2dt_av=(du2dt_p+du2dt_c)*0.50;
du3dt_av=(du3dt_p+du3dt_c)*0.50;
for k=2:n-1
u1(k)=u1_old(k)+du1dt_av(k)*dt;
u2(k)=u2_old(k)+du2dt_av(k)*dt;
u3(k)=u3_old(k)+du3dt_av(k)*dt;
end
%Inlet
u1(1)=rho(1)*A(1);
u2(1)=2*u2(2)-u2(3);
v(1)=(u2(1)/u1(1));
T(1)=(g-1)*(u3(1)/u1(1)-((g/2)*v(1)^2));
u3(1)=u1(1)*(T(1)/(g-1)+(g/2)*v(1)^2);
%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);
m= rho.*A.*v;
p=rho.*T;
rho=u1./A;
v=(u2./u1);
T=(g-1)*((u3./u1)-((g/2)*((u2./u1).^2)));
%Throat Values
rhot2(j)=rho(trt);
vt2(j)=v(trt);
Tt2(j)=T(trt);
%Graph
figure(2)
if j==100
plot(x,m,\'r\',\'linewidth\',1.5);
hold on
elseif j==200
plot(x,m,\'b\',\'linewidth\',1.5);
hold on
elseif j==300
plot(x,m,\'g\',\'linewidth\',1.5);
hold on
elseif j==400
plot(x,m,\'y\',\'linewidth\',1.5);
hold on
elseif j==500
plot(x,m,\'k\',\'linewidth\',1.5);
hold on
elseif j==600
plot(x,m,\'c\',\'linewidth\',1.5);
hold on
elseif j==700
plot(x,m,\'m\',\'linewidth\',1.5);
hold on
elseif j==800
plot(x,m,\'*r\',\'linewidth\',1.5);
hold on
elseif j==900
plot(x,m,\'*b\',\'linewidth\',1.5);
hold on
elseif j==1000
plot(x,m,\'*g\',\'linewidth\',1.5);
grid on
xlabel(\'X\');
ylabel(\'Mass flow rate\');
title(\'Conversative form mass flow rate along with x \',\'color\',\'k\');
legend({\'100 time steps\',\'200 time step\',\'300 time step\',\'400 time step\',\'500 time step\',\'600 time step\',\'800 time step\',\'900 time step\',\'1000 time step\'});
end
end
endc) Main
clear all
close all
clc
%inputs
n=31;
x=linspace(0,3,n);
dx=x(2)-x(1);
%Throat
trt= ((n-1)/2)+1;
%area
A= 1+2.2*(x-1.5).^2
C=0.5
%Gamma
g=1.4;
%Time step
nt=1000;
tic
[rho2,v2,T2,m2,p2,rhot2,vt2,Tt2]= conserve(n,x,dx,A,g,nt,trt,C);
c=toc;
tic
[rho1,v1,T1,m1,p1,rhot1,vt1,Tt1]= non_conver(n,x,dx,A,g,nt,trt,C);
nc=toc;
fprintf(\' Non- conservative time =%d \\n Conservative Time= %d\',nc,c);
%Throat Values
fprintf(\'\\nThroat values for Conservative and Non-Conservative\');
Topic={\'For n\';\'massflowrate\';\'Temperature\';\'Pressure\'};
Convervative=[n;m2(trt);T2(trt);p2(trt)];
NonConvervative= [n;m1(trt);T1(trt);p1(trt)];
table(Topic,Convervative,NonConvervative)
%Exact analytical Solution
m3=linspace(0.579,0.579,n);
z=linspace(1,nt,nt);
figure(3)
plot(x,m1,\'r\',\'linewidth\',1.5);
hold on
plot(x,m2,\'b\',\'linewidth\',1.5);
hold on
plot(x,m3,\'- k\');
grid on
xlabel(\'X\');
ylabel(\'Mass flow rate\');
title(\'Mass flow comparision between Conservative and Non-conservative \',\'color\',\'k\');
legend({\'Non-conservative\',\'Conservative\',\'Exact Analytical\'});
%Throat Values
figure(4)
subplot(3,1,1);
plot(z,rhot2,\'b\');
hold on
plot(z,rhot1,\'r\');
grid on
xlabel(\'timesteps\');
ylabel(\'rho\');
title(\'Throat Values comparision between Conservative and Non-conservative (n=31)\',\'color\',\'k\');
legend({\'Conservative\',\'Non-conservative\'});
subplot(3,1,2);
plot(z,vt2,\'b\');
hold on
plot(z,vt1,\'r\');
grid on
xlabel(\'timesteps\');
ylabel(\'V\');
legend({\'Conservative\',\'Non-conservative\'});
subplot(3,1,3);
plot(z,Tt2,\'b\');
hold on
plot(z,Tt1,\'r\');
grid on
xlabel(\'timesteps\');
ylabel(\'T\');
legend({\'Conservative\',\'Non-conservative\'});Results:
For n=31
1) Non-conservative mass flow vs distance through the nozzle
It can be seen that after 600-time steps, the mass flow distribution is a straight, horizontal line across the graph. Therefore, the mass flow rate has converged to a constant steady-state value throughout the nozzle.
2) Conservative mass flow vs distance through the nozzle
It can be seen that after 800-time steps, the mass flow distribution is a straight, horizontal line across the graph. Therefore, the mass flow rate has converged to a constant steady-state value throughout the nozzle.
3) steady-state mass flow variations comparison between conservative and non-conservative forms
This figure shows that conservative form performs better mass preservation throughout the nozzle. Its graph is more closer to the Exact analytical graph. However, at the throat, Non-conservative gives more accurate results.
4)Timewise Variation of certain factors at nozzle throat:
For n=61
5) Timewise Variation of certain factors at nozzle throat:
It can be seen from the graph that the graphs are similar only the time steps taken for each activity have doubled. This is because we have doubled the grid points. It can be viewed as the two steady-state solutions that are independent of grid points. It can be further explained by the following figures. Also, it can be seen that for n=61, the solution does not converge at 1000 time steps, therefore, calculating n=61 at 1600 time steps we get;
Therefore, it can be seen that as the no. of the grid points are increasing the values are coming close to that of analytical results( rho=0.634). Hence, the numerical solution is improving but it is doing so only marginally for both the forms, which indicates grid independence.
Also, it should be noted that the results of the non-conservation form are closer to the Analytical solution that that of conservation form.
Faster Method:
It can be seen that the non-conservative method is faster as it does not have any flux variables calculation as in the conservative method.
[References from Computational Fluid Dynamics by John D. Anderson be taken for writing theory]