1D Super-sonic nozzle flow simulation using Macormack Method
AIM
- Write code for the 1D supersonic nozzle flow equations using the Macormack Method, for conservative and non-conservative approach.
- Figure out the minimum number of cycles for which the simulation should be run in order for convergence
- Compare the normalized mass flow rate between the conservative forms and non-conservative forms - make suitable conclusions
Introduction
For 1D supersonic flow through the convergent-divergent nozzle, steady-state isentropic flow is considered. Inlet of the nozzle is subsonic and at the throat is sonic and at the outlet is supersonic.
The length of the nozzle is taken as 3 units, and it is firstly discretized using 31 units. Tolerance for convergence for both conservative and non-conservative is taken as 1e-4.
The solution is solved using the time-marching solution using the MacCormack method.
Maccormack Method
It is explicit finite difference techniques which are second-order accurate in both space and time. It is very effective where we can obtain second-order accuracy with the first-order terms and mathematical complexity of calculating the second-order time, derivatives can be avoided. In CFD, the Maccromack method is a widely used discretization scheme for numerical calculation of hyperbolic partial differential equations.
According to the MacCormack method, the value of a variable at grid point i,j at time level t+dt is given as follows...
`v_(i,j)^(t+trianglet) = v_(i,j)^(t) + (frac{delv}{delt})_(avg).trianglet`
where the average time derivatives are the mean value between time t and t+dt. The average time derivatives are obtained from the predictor-corrector philosophy.
`(frac{delv}{delt})_(avg) = 0.5[(frac{delv}{delt})_(predicted) + (frac{delv}{delt})_text(corrected)]`
Predicted Step
in the predictor step, the spatial term on the right-hand side of the equation is replaced by forwarding difference. The RHS of the equation is known as all the value of the variable are for time-step t. Next, the value of the variable for LHS is calculated for time step t+dt.
Consider the following equations,
`frac{delv}{delt} = frac{delrho}{delx}`
`(frac{delv}{delt})_text(predicted) = (frac{rho_text(i+1,j)^t - rho_text(i,j)^t }{trianglex} +frac{rho_text(i+1,j)^t - rho_text(i,j)^t }{triangley} )`
`v_text(i,j)^(t+trianglet)=v_text(i,j)^(t)+(frac{delv}{delt})_text(predicted). trianglet`
Corrected Step
In the corrected step, the spatial term on the right side are replaced by a backward difference. The value of the variable is the RHS are taken which are predicted step for time step 't+dt'. The calculation for corrected time derivatives is.
`(frac{delv}{delt})_text(corrected) = (frac{rho_text(i,j)^(t+trianglet) - rho_text(i-1,j)^(t+trianglet) }{trianglex} +frac{rho_text(i,j)^(t+trianglet) - rho_text(i-1,j)^(t+trianglet) }{triangley} )`
Now, the corrected and predicted value of time derivatives is used to calculate the average time derivatives which in turn gives the value of the variable at a given point for a particular time-step, it's repeated for all grid points.
Governing Equations
Conservative Form
The conservative form of expression can be expressed as consists of vector U, flux vector F, and source term J.
`frac{delU}{delt} + frac{delF}{delx} = J`
the equation of quasi 1-D supersonic flow can be expressed in a similar generic form as follows
`U_1 = rhoAF_1 = rhoAu`
`U_2= rhoAu`
`F_2 = rhoAu^2 +frac{1}{gamma}pA`
`J_2 = frac{1}{gamma}pfrac{delA}{delx}`
`U_3 = rho(frac{e}{gamma-1} + frac{gamma}{2}u^2)`
`F_3 = rho(frac{e}{gamma-1} + frac{gamma}{2}u^2) + frac{1}{gamma}pA`
The equation which is solved using this Maccromack method are-
`frac{delU_1}{delt} = frac{delF_1}{delt}`
`frac{delU_2}{delt} = -frac{delF_2}{delx} + J_2`
`frac{delU_3}{delt} = -frac{delF_3}{delx}`
further, the flux vector and source term are expressed in terms of solution vectors. Also to find the field variables, they are written in terms of solution vectors.
`F_1 = U_1`
`F_2 = frac{U_2^2}{U_1} + frac{gamma-1}{gamma}(U_3 - frac{gamma}{2}frac{U_2^2}{U_1})`
`F_3 = gammafrac{U_2U_3}{U_1} + frac{gamma(gamma-1)}{2}frac{U_2^3}{U_1^2}`
`J_2 = frac{gamma(gamma-1)}{gamma}(U_3 - frac{gamma}{2}frac{U_2^2}{U_1})frac{del lnA}{delx}`
Also,
`rho = frac{U_1}{A} , u = frac{U_2}{U_1} , T = (gamma-1)(frac{U_3}{U_1} - frac{gamma}{2}u^2)`
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);
Initial Conditions
Area = 1+(2.2*(x-1.5).^2);
Velocity= 0.59./(rho.*A);
`rho = 1, T = 1; text(for) ,0<=x<=0.5`
`rho = 1-(0.366*(x-0.5)); T = 1-(0.167*(x-0.5)); x >=0.5 & x < 1.5`
`rho = 0.634-(0.3879*(x-1.5)); T = 0.833-(0.3507*(x-1.5)); x>=1.5 & x < 3.5`
Non-Conservative form
Continuity Equation
`frac{delrho}{delt} = -rhofrac{delu}{delx} - rhoufrac{del(lnA)}{delx} - ufrac{delrho}{delx}`
Momentum Equation
`frac{delu}{delt} = -ufrac{delu}{delx} - frac{1}{gamma}(frac{delT}{delx} + frac{T}{rho}frac{delrho}{delx})`
Energy Equation
`frac{delT}{delt} = -ufrac{delT}{delx} -(gamma-1)T(frac{delu}{delx} + ufrac{del(lnA)}{delx})`
Calculation of time-step
For a linear hyperbolic equation, the value of CFL<1, will give a stable solution. And the time-step is calculated as:
`Deltat = C frac{trianglex}{V+a}`
where, C is Courant Number
V is local flow velocity at a point
a is the local speed of sound
Boundary Conditions
Inlet
Velocity at the inlet is given by,
`u_1 = 2u_2 + u_3`
Outlet
At the outlet, variable are keep changing with respect to time, and it is given by
`u_n = 2u_(n-1) + u_(n-2)`
`T_n = 2T_(n-1) + T_(n-2)`
`rho_n = 2rho_(n-1) + rho_(n-2)`
Initial Conditions
Area of nozzle will remain constant throughout the time of simulation
Area = 1+(2.2*(x-1.5).^2);
other flow variables are given by,
`rho = 1+ 0.3146x`
`T = 1+ 0.2314x`
`u = (0.1 + 1.09x)T^0.5`
Matlab Code
Main Method
clc
clear all
close all
n = 31;
x = linspace(0,3,n);
dx = x(2)-x(1);
gamma = 1.4;
C = 0.5;% courant number
throat = ((n-1)/2) + 1;
% non-dimensional area
A = 1+(2.2*(x-1.5).^2);
% time-step
nt = 1400;
dt = 0.01;
%non-conservative-form
tic
[rho,T,v,P,M,m,rho_throat,T_throat,v_throat,P_throat,M_throat,m_throat] = non_conservative_form(n,x,dx,gamma,throat,nt,dt,A);
time1 = toc;
%conservative-form
tic
[rho_c,T_c,v_c,P_c,M_c,m_c,rho_throat_c,T_throat_c,v_throat_c,P_throat_c,M_throat_c,m_throat_c] = conservative_form(n,x,dx,gamma,throat,nt,dt,A);
time2=toc;
figure(3)
subplot(2,3,1)
plot(x,rho,'r',x,rho_c,'b*');
xlabel('Length of Nozzle,m');
ylabel('Non-dimensional density');
title('Density-Profile');
legend('Non-Conservative','Conservative');
subplot(2,3,2)
plot(x,T,'r',x,T_c,'b*');
xlabel('Length of Nozzle,m');
ylabel('Non-Dimensinal Temperature');
title('Non-Conservative,Conservative');
subplot(2,3,3)
plot(x,v,'r',x,v_c,'b*');
xlabel('Length of Nozzle,m');
ylabel('Non-dimensional Velocity');
title('Velocity-Profile');
legend('Non-Conservative','Conservative');
subplot(2,3,4)
plot(x,P,'r',x,P_c,'b*');
xlabel('Length of Nozzle,m');
ylabel('Non-dimensional Pressure');
title('Pressure-Profile');
legend('Non-Conservative','Conservative');
subplot(2,3,5)
plot(x,M,'r',x,M_c,'b*');
xlabel('Length of Nozzle,m');
ylabel('Mach Number');
title('Mach Number Profile');
legend('Non-Conservative','Conservative');
subplot(2,3,6)
plot(x,m,'r',x,m_c,'b*');
xlabel('Length of Nozzle,m');
ylabel('Mass Flow Rate');
title('Mass Flow Rate Profile');
legend('Non-Conservative','Conservative');
figure(4)
subplot(2,3,1)
plot(1:nt,rho_throat,'r',1:nt,rho_throat_c,'k');
hold on
line([0 nt],[0.634 0.634],'color','b');
xlabel('time,s');
ylabel('Non-dimensiona Density');
title('Variation of density at Nozzle throat')
grid minor
legend('Non-Conservative','Conservative','Exact Solution');
subplot(2,3,2)
plot(1:nt,T_throat,'r',1:nt,T_throat_c,'k');
hold on
line([0 nt],[0.833 0.833],'color','b');
xlabel('time,s');
ylabel('Non-dimensiona temperature');
title('Variation of temperature at Nozzle throat')
grid minor
legend('Non-Conservative','Conservative','Exact Solution');
subplot(2,3,3)
plot(1:nt,v_throat,'r',1:nt,v_throat_c,'k');
xlabel('time,s');
ylabel('Non-dimensional Velocity');
title('Variation of Velocity at Nozzle throat')
grid minor
legend('Non-Conservative','Conservative');
subplot(2,3,4)
plot(1:nt,P_throat,'r',1:nt,P_throat_c,'k');
hold on
line([0 nt],[0.528 0.528],'color','b');
xlabel('time,s');
ylabel('Non-dimensional Pressure');
title('Variation of Pressure at Nozzle throat')
grid minor
legend('Non-Conservative','Conservative','Exact Solution');
subplot(2,3,5)
plot(1:nt,M_throat,'r',1:nt,M_throat_c,'k');
hold on
line([0 nt],[1 1],'color','b');
xlabel('time,s');
ylabel('Mach Number');
title('Variation of Mach Number at Nozzle throat')
grid minor
legend('Non-Conservative','Conservative','Exact Solution');
subplot(2,3,6)
plot(1:nt,m_throat,'r',1:nt,m_throat_c,'k');
hold on
line([0 nt],[0.579 0.579],'color','b');
xlabel('time,s');
ylabel('Non-dimensional Mass Flow Rate');
title('Variation of Mass Flow Rate at Nozzle throat')
grid minor
legend('Non-Conservative','Conservative','Exact Solution');
Conservative Method
function [rho,T,v,P,M,m,rho_throat,T_throat,v_throat,P_throat,M_throat,m_throat] = conservative_form(n,x,dx,gamma,throat,nt,C,A)
% initial profile
for i = 1:n
if x(i)>=0 && x(i)< 0.5
rho(i) = 1;
T(i) = 1;
else
if 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));
else
if 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
end
end
% velocity and pressure
v = 0.59./(rho.*A);
P = rho.*T;
% U1 ,U2 and U3
U1 = rho.*A;
U2 = rho.*A.*v;
U3 = rho.*A.*((T/(gamma-1))+(0.5*gamma.*v.*v));
% time loop
for k = 1:nt
U1_init = U1;
U2_init = U2;
U3_init = U3;
% time-step
dt = min(0.5.*dx./(sqrt(T) + v));
%flux vector
F1 = U2;
F2 = ((U2.^2)./(U1))+((gamma-1)/gamma).*(U3-(0.5*gamma.*U2.^2)./U1);
F3 = ((gamma.*U2.*U3)./U1) - (0.5*gamma*(gamma-1).*U2.^3)./(U1.^2);
%predicted steps
for i = 2:n-1
dF1_dx = (F1(i+1)-F1(i))/dx;
dF2_dx = (F2(i+1)-F2(i))/dx;
dF3_dx = (F3(i+1)-F3(i))/dx;
%continuity equations
dU1_dt_p(i) = -dF1_dx ;
%momentum equations
J2(i) = (1/gamma).*rho(i).*T(i).*((A(i+1)-A(i))/dx);
dU2_dt_p(i) = -dF2_dx + J2(i) ;
% energy equations
dU3_dt_p(i) = -dF3_dx ;
%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
% update
rho = U1./A;
v = U2./U1;
T = (gamma-1).*((U3./U1)-((gamma*v.^2)/2));
P = rho.*T;
F1 = U2;
F2 = ((U2.^2)./(U1))+((gamma-1)/gamma).*(U3-(0.5*gamma.*U2.^2)./U1);
F3 = ((gamma.*U2.*U3)./U1) - (0.5.*gamma.*(gamma-1).*U2.^3)./(U1.^2);
% correction-step
for j = 2:n-1
dF1_dx = (F1(j)-F1(j-1))/dx;
dF2_dx = (F2(j)-F2(j-1))/dx;
dF3_dx = (F3(j)-F3(j-1))/dx;
% continuity equations
dU1_dt_c(j) = -dF1_dx ;
%momentum equations
J2(j) = (1/gamma).*rho(j).*T(j).*((A(j)-A(j-1))/dx);
dU2_dt_c(j) = -dF2_dx + J2(j) ;
% energy equations
dU3_dt_c(j) = -dF3_dx ;
end
% average time derivatives
dU1_dt_avg = 0.5.*(dU1_dt_p +dU1_dt_c);
dU2_dt_avg = 0.5.*(dU2_dt_p +dU2_dt_c);
dU3_dt_avg = 0.5.*(dU3_dt_p +dU3_dt_c);
% final solution update
for l = 2:n-1
U1(l) = U1_init(l) +dU1_dt_avg(l)*dt;
U2(l) = U2_init(l) +dU2_dt_avg(l)*dt;
U3(l) = U3_init(l) +dU3_dt_avg(l)*dt;
end
% apply 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);
%update
rho = U1./A;
v = U2./U1;
T = (gamma-1).*((U3./U1)-((gamma*v.^2)/2));
P = rho.*T;
%mach number and mass flow rate
M = v./((T).^0.5);
m = rho.*v.*A;
% value sat throat
rho_throat(k) = rho(throat);
v_throat(k) = v(throat);
T_throat(k) = T(throat);
P_throat(k) = P(throat);
M_throat(k) = M(throat);
m_throat(k) = m(throat);
%plot
b = figure(2);
hold on
if k ==100
plot(x,m,'b','linewidth',2);
else
if k ==150
plot(x,m,'g','linewidth',2);
else
if k ==200
plot(x,m,'m','linewidth',2);
else
if k ==700
plot(x,m,'r','linewidth',2);
legend('100delta','150delta','200delta','700delta');
title('Variation of Mass flow Rate at different time-step');
xlabel('Non-dimensional length of Nozzle')
ylabel('Non-dimensional mass flow rate')
grid minor
saveas(b,'con_mass.png')
end
end
end
end
endNon-Conservative
function [rho,T,v,P,M,m,rho_throat,T_throat,v_throat,P_throat,M_throat,m_throat] = non_conservative_form(n,x,dx,gamma,throat,nt,C,A)
% Calculating initial profile
rho = 1-(0.3146*x);
T = 1-(0.2314*x);
v = (0.1+(1.09*x)).*T.^(0.5);
%outer-time-loop
for k = 1:nt
rho_ini = rho;
T_ini = T;
v_ini = v;
%time-step
dt = 0.01;
%predicted -step
for i = 2:n-1
dv_dx = (v(i+1)-v(i))/dx;
drho_dx = (rho(i+1)-rho(i))/dx;
dlogA_dx = (log(A(i+1))-log(A(i)))/dx;
dT_dx = (T(i+1)-T(i))/dx;
%continuity-equation
drho_dt_p(i) = -rho(i)*dv_dx-rho(i)*v(i)*dlogA_dx-v(i)*drho_dx;
%momentum equation
dv_dt_p(i) = -v(i)*dv_dx-(1/gamma)*(dT_dx+((T(i)/rho(i))*drho_dx));
%energy equation
dT_dt_p(i) = -v(i)*dT_dx-(gamma-1)*T(i)*(dv_dx+(v(i)*dlogA_dx));
%update
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
% corrected step
for j = 2:n-1
dv_dx_c = (v(j)-v(j-1))/dx;
drho_dx_c = (rho(j)-rho(j-1))/dx;
dlogA_dx_c = (log(A(j))-log(A(j-1)))/dx;
dT_dx_c = (T(j)-T(j-1))/dx;
%continuity-equation
drho_dt_c(j) = -rho(j)*dv_dx_c-rho(j)*v(j)*dlogA_dx_c-v(j)*drho_dx_c;
%momentum equation
dv_dt_c(j) = -v(j)*dv_dx_c-(1/gamma)*(dT_dx_c+((T(j)/rho(j))*drho_dx_c));
%energy equation
dT_dt_c(j) = -v(j)*dT_dx_c-(gamma-1)*T(j)*(dv_dx_c+(v(j)*dlogA_dx_c));
end
%average time derivatives
drho_dt = 0.5.*(drho_dt_p + drho_dt_c);
dv_dt = 0.5.*(dv_dt_p + dv_dt_c);
dT_dt = 0.5.*(dT_dt_p + dT_dt_c);
%solution update
for l = 2:n-1
rho(l) = rho_ini(l) + drho_dt(l)*dt;
T(l) = T_ini(l) +dT_dt(l)*dt;
v(l) = v_ini(l) + dv_dt(l)*dt;
end
%boundary conditions
%inlet
v(1) = 2*v(2)-v(3);
rho(1) = 1;
T(1) = 1;
%outlet
rho(n) = 2*rho(n-1) - rho(n-2);
v(n) = 2*v(n-1) - v(n-2);
T(n) = 2*T(n-1) - T(n-2);
% pressure,machnumber,mass-flow-rate
P = rho.*T;
M = (v./(T).^0.5);
m = rho.*v.*A;
%throat values
rho_throat(k) = rho(throat);
v_throat(k) = v(throat);
T_throat(k) = T(throat);
P_throat(k) = P(throat);
M_throat(k) = M(throat);
m_throat(k) = m(throat);
%plotting
a = figure(1);
hold on
if k ==100
plot(x,m,'b','linewidth',2);
else
if k ==150
plot(x,m,'g','linewidth',2);
else
if k ==200
plot(x,m,'m','linewidth',2);
else
if k ==700
plot(x,m,'r','linewidth',2);
legend('100delta','150delta','200delta','700delta');
title('Variation of Mass flow Rate at different time-step');
xlabel('Non-dimensional length of Nozzle')
ylabel('Non-dimensional mass flow rate')
grid minor
saveas(a,'non-con_mass.png')
end
end
end
end
end
endRESULT
1. Variation of FLow along nozzle Length
From the above plot, both the conservative and non-conservative approach is in close proximity to every parameter. And we can also see that temperature, density, and pressure decreases along the length, whereas velocity and Mach number decreases.
2. Variation of Flow along the throat
From the above plot, we can conclude that conservative and non-conservative have become stable at around 700 time-step.
We also draw a conclusion, from the plot that the non-conservative approach is in close proximity to the analytical solution.
Conservative solution approaches stability earlier than a non-conservative approach.
3. The mass-flow rate at different time-step
Conservative Approach
From the above graph, we can see that as the number of time-step increases the solution becomes more stable and fluctuation decreases.
As we can see fro time-step 600, and 700 they are stable and are in close proximity.
Non-Conservative Approach
From the above graph, we can see that as the number of time-step increases the solution becomes more stable and fluctuation decreases.
As we can see fro time-step 600, and 700 they are stable and are in close proximity.
Grid Dependency Test
For the grid dependency test, the solution was tested at two different grid for 31 and 61. Result we obtain is shown below
From the above solution, we can conclude that the solution is independent of the grid point. The solution almost remains the same at both grid points. We can't see any change in pressure, temperature, and density at both grid points.
Conclusion
For grid point 31, and time step 1400, the computation time required for a different approach is:
The computational time for the non-conservative approach is: 8.0433 seconds
The computational time for the conservative approach is: 5.2312 seconds
Hence from this finding, we can conclude that the conservative approach is faster than the non-conservative approach.