Week 7 - Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
The objective of the Problem :
Simulate isentropic flow through a quasi 1D subsonic-supersonic nozzle using both the conservation and non-conservation forms of the governing equations, solve them using MacCormack's approach, and compare the results.
Explanation :
We consider the study, isentropic flow through a convergent-divergent nozzle. The flow from the inlet to the nozzle comes from
the reservoir. The cross-section area of the reservoir is very large and hence velocity is very small. The flow is locally subsonic in the convergent
section of the nozzle, Sonic at the throat and Supersonic at the divergent section.
Consider a small control volume :
We assume that at a given section, where the cross-section area is A, the flow properties are uniform across the section. Hence, although the
area of the nozzle changes as a function of distance along with the nozzle, X, and therefore in reality the flow field is two dimensional(the flow
varies in the two-dimensional x-y space), We make the assumption that the flow properties vary only with X; this is tantamount to assuming
uniform flow properties across any given cross-section. Such flow is defined as Quasi one dimensional flow.
The governing continuity, momentum and energy equations in terms of non-dimensional variables for this quasi-one-dimensional steady
and isentropic flow can be expressed respectively as:
Non Conservative continuty equation :
`frac(delrho')(delt') = -rho'(frac(delV')(delx') - rho'V'(frac(del(loga'))(delx'))-V'(frac(delrho')(delx'))`
Non Conservative momentum equations :
`frac(delV')(delt') = -V'(frac(delV')(delx')) - frac(1)(gamma)(frac(delT')(delx')+frac(T')(rho')(frac(delrho')(delx')))`
Non Conservatice energy equation :
`frac(delT')(delt') = -V'(frac(delT')(delx') )- (gamma-1)T'[frac(delV')(delx')+V'frac(del(logA'))(delx')]`
Conservative continuty Equation :
`frac(rho'A')(delt') + frac(del(rho'A'V'))(delx') = 0`
Conservative Momemtum equation :
`frac(del(rho'A'V'))(delt') + frac(del[rho'A'(V')^2+(frac(1)(gamma))rho'A'])(delx') = frac(1)(gamma)rho'(frac(delA')(delx'))`
Conservation energy equation :
`frac(del[rho'(frac(T')(gamma-1)+frac(gamma)(2)V^2)A'])(delt') + frac(del[rho'(frac(T')(gamma-1)+frac(gamma)(2)V^2)V'A'+rho'A'V'])(delx') = 0`
McCormack method:
To discretize the numerical solution we use the technique called MacCormack method. The MacCormack method is a variation of the two- step Lax-Wendroff Scheme but is much simpler in the application.
The application of the MacCormack method proceeds in two steps; a predictor step which is followed by a corrector step.
Initial Profile :
Before the calculation the initial profile for all the flow field variables which are velocity, density and temperature have to be determined. These are guess values but should be specified with more awareness on the problem
Predictor step: In the predictor step, a "provisional" value of u at timelevel n + 1 (denoted by `(u_i)^(n+1)`) is estimated obtained by replacing thespatial and temporal derivatives in the previous first-order equation using Forward Difference Scheme. `Deltat` is used for Timestep.
`(frac(delf)(delt))_p = frac(f(i+1)-f(i))(Deltax)`
solution update :
`f_p = f_(Initial) + (frac(delf)(delt))_p.Deltat`
`f_p` = Flow feild variables calculated after yhr predictor step
`f_(Inital)` = flow feild variable in the initial step
`Deltat`time step size
Corrector step: In the corrector step, the predicted value `(u_i)^(n+1)` corrected according to the equation. Note that the corrector step uses the Backward Difference Scheme for spatial derivative. The time-step used in the corrector step is `frac(Deltat)(2)` in contrast to the `Deltat` used in the predictor step.
`(frac(delf_p)(delt))_c = frac(f(i) - f(i-1))(Deltax)`
`(frac(delf)(delt))_(avg) = 0.5((frac(delf)(delt))_p+(frac(delf_p)(delt))_c)`
Final Solution update :
`f_(sol) = f_(Initial) + (frac(delf)(delt))_(avg).Deltat`
where `f_(sol)`= Final solution of the flow feild variables for the first time step
Main Code :
clear all
close all
clc
%Input Variables
L=3;%length of nozzle
n =31;%Number of grid points
x=linspace(0,L,n);%mesh definition
dx = x(2)-x(1);%element size
throat = (n+1)/2 %throat position
nt=1400; %number of iterations
a = 1+2.2.*(x-1.5).^2; %area
C= 0.5; %CFL Number
gamma = 1.4; %cp/cv
%function of Non conservative method
[rho_NC,v_NC,T_NC,P_NC,M_NC,Mass_flow_NC]=Non_Conservative_function(n,x,dx,gamma,throat,nt,a,C);
%function of conservative method
[rho_C,v_C,T_C,P_C,M_C,Mass_flow_C]=Conservative_function(n,x,dx,gamma,throat,nt,a,C);
Terms = [{'density'};{'Tempearture'};{'pressure'};{'Mach number'};{'Velocity'}];
Exact_Sol = [0.634;0.833;0.528;1;1];
Non_Conservative_Sol = [rho_NC(throat);T_NC(throat);P_NC(throat);M_NC(throat);v_NC(throat)];
Conservative_Sol = [rho_C(throat);T_C(throat);P_C(throat);M_C(throat);v_C(throat)];
table(Terms,Exact_Sol,Non_Conservative_Sol,Conservative_Sol)
figure(7)
plot(x,Mass_flow_NC,'m','linewidth',1)
hold on
plot(x,Mass_flow_C,'y','linewidth',1)
xlabel('length of nozzel')
ylabel('mass flow rate')
legend('Non-conservative form','Conservative form')
title('Comparison mass flow distributions of Non conservatiove form & Conservative form')
grid onExplanation:
• Code is started with clear all, close all and clc to enable a fresh start without the manipulation of previous values.
• Input parameters are assigned for the variables x, dx, n, gamma, nt and C.
• Non-Conservative and Conservative forms are created as separate functions with the necessary inputs which is detailed below this
section.
• tic to command is used for calculating the time elapsed between the functions.
• Plotting is done for the following solutions
1) Timestep variation at throat area of both forms
2) Variation of Flow Rate Distribution of both forms
3) Comparison of Normalized Mass Flow Rate of both forms
• Each plot is computed as subplots of four primitive variables such as Density, Temperature, Pressure and Mach Number at the throat as
well as Convergent/Divergent section.
Function Code: Non-conservative form
function [rho_NC,v_NC,T_NC,P_NC,M_NC,Mass_flow_NC]=Non_Conservative_function(n,x,dx,gamma,throat,nt,a,C);
rho_NC = 1-0.31468*x; %Density
T_NC = 1-0.2341*x; %Temperature
v_NC = (0.1+1.09*x.*T_NC.^0.5); %Velocity
tic
for k = 1:nt
rho_old = rho_NC;
T_old = T_NC;
v_old = v_NC;
dt = min(C.*dx./(sqrt(T_NC)+v_NC));
for j = 2:n-1
dv_dx = (v_NC(j+1) - v_NC(j))/dx;
drho_dx = (rho_NC(j+1) - rho_NC(j))/dx;
dT_dx = (T_NC(j+1) - T_NC(j))/dx;
dloga_dx = (log(a(j+1)) - log(a(j)))/dx;
drho_dt_P(j) = -rho_NC(j)*dv_dx - rho_NC(j)*v_NC(j)*dloga_dx - v_NC(j)*drho_dx;
dv_dt_P(j) = -v_NC(j)*dv_dx - (1/gamma)*(dT_dx + (T_NC(j)*drho_dx/rho_NC(j)));
dT_dt_P(j) = -v_NC(j)*dT_dx - (gamma-1)*T_NC(j)*(dv_dx + v_NC(j)*dloga_dx);
v_NC(j) = v_NC(j) + dv_dt_P(j)*dt;
rho_NC(j) = rho_NC(j) + drho_dt_P(j)*dt;
T_NC(j) = T_NC(j) + dT_dt_P(j)*dt;
end
for j = 2:n-1
dv_dx = (v_NC(j) - v_NC(j-1))/dx;
drho_dx = (rho_NC(j) - rho_NC(j-1))/dx;
dT_dx = (T_NC(j) - T_NC(j-1))/dx;
dloga_dx = (log(a(j)) - log(a(j-1)))/dx;
drho_dt_C(j) = -rho_NC(j)*dv_dx - rho_NC(j)*v_NC(j)*dloga_dx - v_NC(j)*drho_dx;
dv_dt_C(j) = -v_NC(j)*dv_dx - (1/gamma)*(dT_dx + (T_NC(j)*drho_dx/rho_NC(j)));
dT_dt_C(j) = -v_NC(j)*dT_dx - (gamma-1)*T_NC(j)*(dv_dx + v_NC(j)*dloga_dx);
end
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);
for i = 2:n-1
rho_NC(i) = rho_old(i) + drho_dt(i)*dt;
v_NC(i) = v_old(i) + dv_dt(i)*dt;
T_NC(i) = T_old(i) + dT_dt(i)*dt;
end
v_NC(1) = 2*v_NC(2) - v_NC(3);
v_NC(n) - 2*v_NC(n-1)-v_NC(n-2);
rho_NC(n) = 2*rho_NC(n-1) - rho_NC(n-2);
T_NC(n) = 2*T_NC(n-1) - T_NC(n-2);
P_NC = rho_NC.*T_NC;
M_NC = (v_NC./sqrt(T_NC));
Mass_flow_NC = rho_NC.*v_NC.*a;
rho_throat(k) = rho_NC(throat);
v_throat(k) = v_NC(throat);
T_throat(k) = T_NC(throat);
P_throat(k) = P_NC(throat);
M_throat(k) = M_NC(throat);
flow_rate_throat(k) = Mass_flow_NC(throat);
figure(1);
if k==100
plot(x,Mass_flow_NC,'m','linewidth',1)
hold on
else if k == 200
plot(x,Mass_flow_NC,'y','linewidth',1)
hold on
elseif k ==500
plot(x,Mass_flow_NC,'k','linewidth',1)
hold on
else if k == 1000
plot(x,Mass_flow_NC,'w','linewidth',1)
else if k == 1200
plot(x,Mass_flow_NC,'b','linewidth',1)
hold on
legend('100 Delta','200 Delta','500 Delta','1000 Delta','1200 Delta')
title('Mass flow rate distribution at different time steps in Non-conservation form')
xlabel('Non dimensional length of nozzle')
ylabel('Mass flow rate in Non-Conservative from')
grid on
end
end
end
end
end
Time_NC = toc
figure(2)
plot(M_throat,'m','linewidth',1)
hold on
plot(rho_throat,'y','linewidth',1)
hold on
plot(v_throat,'k','linewidth',1)
hold on
plot(T_throat,'w','linewidth',1)
hold on
plot(P_throat,'b','linewidth',1)
legend('Mach number','(rho/rho_o)','(v/a_o)','(T/T_o)','(P/P_o)')
xlabel('No.of Time steps')
ylabel('M,rho/rho_o,v/v_o,T/T_o,P/P_o')
title('Timewise variation non dimensional flow feild variables at throat in non-Conservative form')
grid on
figure(3)
subplot(5,1,1)
plot(x,M_NC,'m','linewidth',1)
ylabel({'mach number'})
xlabel({'macch number at throat'})
subplot(5,1,2)
plot(x,v_NC,'y','linewidth',1)
ylabel({'Velocity'})
xlabel({'Velocity at throat'})
subplot(5,1,3)
plot(x,rho_NC,'k','linewidth',1)
ylabel({'Density'})
xlabel({'Density at throat'})
subplot(5,1,4)
plot(x,T_NC,'w','linewidth',1)
ylabel({'Temperature'})
xlabel({'Temperature at throat'})
subplot(5,1,5)
plot(x,P_NC,'b','linewidth',1)
ylabel({'pressure'})
xlabel({'pressure at throat'})
endExplanation:
The purpose of the function is to calculate the values of
Mass flow rate, Pressure, Mach number, Density, velocity, Temperature, Density at the throat, the velocity at the throat, Temperature at the throat,
Mass flow rate at the throat, Pressure at the throat, Mach number at the throat by the known equations for all the grid points.
• For this initial boundary conditions are established for Density, Temperature & Velocity and profiles for calculating area & throat.
• Now Outer Time loop is computed for loop to calculate all the grid points for various formulas inside the loop.
• MacCarmack technique is used to solve the numerical solution of governing equations such as Continuity, Momentum & Energy
Equation.
The technique implies two methods for solving namely the Predictor method followed by the Corrector method.
• Predictor & Corrector method are assigned within the for loop individually for calculating values.
• The calculated values are averaged in the next step and updating the solution with old values.
• Boundary Conditions are applied within the for loop and subsequently used for calculating the Mass Flow Rate, Mach
Number & Pressure and also calculating the values throat section.
• Finally plotting is done for Distance (x) & Mass Flow Rate for different Timesteps.
• Program is ended.
Function Code: Conservative form
function [rho_C,v_C,T_C,P_C,M_C,Mass_flow_C]=Conservative_function(n,x,dx,gamma,throat,nt,a,C);
tic
for i = 1:n
if x(i) <0.5
rho_C(i) = 1;
T_C(i) = 1;
elseif (x(i) >=0.5 && x(i) <=1.5)
rho_C(i)= 1-0.366*(x(i)-0.5);
T_C(i) = 1-0.167*(x(i)-0.5);
elseif (x(i) >=1.5 && 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
v_C(i) = 0.59/(rho_C(i)*a(i));
end
dt = min(C.*dx./(sqrt(T_C)+v_C));
%dt=0.0267;
U1 = rho_C.*a;
U2 = rho_C.*a.*v_C;
U3 = U1.*(T_C/(gamma - 1)) + ((gamma/2).*v_C.^2);
tic
for k = 1:nt
U1_old = U1;
U2_old = U2;
U3_old = U3;
F1 = U2;
F2 = U2.^2./U1 + ((gamma -1)/gamma).*(U3 - (gamma/2)*U2.^2./U1);
F3 = gamma.*U2.*U3./U1 - 0.5*gamma*(gamma-1).*U2.^3./U1.^2;
for i = 2:n-1
J2 = (1/gamma)*rho_C(i)*T_C(i)*((a(i+1)-a(i))/dx);
dU1_dt_P(i) = -((F1(i+1) - F1(i))/dx);
dU2_dt_P(i) = -((F2(i+1) - F2(i))/dx) + J2;
dU3_dt_P(i) = -((F3(i+1) - F3(i))/dx);
U1(i) = U1_old(i) + dU1_dt_P(i)*dt;
U2(i) = U2_old(i) + dU2_dt_P(i)*dt;
U3(i) = U3_old(i) + dU3_dt_P(i)*dt;
end
F1 = U2;
F2 = U2.^2./U1 + ((gamma -1)/gamma).*(U3 - (gamma/2)*U2.^2./U1);
F3 = gamma*U2.*U3./U1 - 0.5*gamma*(gamma-1).*U2.^3./U1.^2;
rho_C = U1./a;
T_C = (gamma-1).*(U3./U1-(gamma/2).*(U2./U1).^2);
for i = 2:n-1
J2 = (1/gamma)*rho_C(i)*T_C(i)*((a(i)-a(i-1))/dx);
dU1_dt_C(i) = -((F1(i) - F1(i-1))/dx);
dU2_dt_C(i) = -((F2(i) - F2(i-1))/dx) + J2;
dU3_dt_C(i) = -((F3(i) - F3(i-1))/dx);
end
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);
for j = 2:n-1
U1(j) = U1_old(j) + (dU1_dt(j))*dt;
U2(j) = U2_old(j) + (dU2_dt(j))*dt;
U3(j) = U3_old(j) + (dU3_dt(j))*dt;
end
U1(1) = a(1);
U2(1) = 2*U2(2) - U2(3);
v_C(1) = U2(1)/U1(1);
U3(1) = U1(1)*((T_C(1)/(gamma-1))+((gamma/2)*(v_C(1))^2));
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);
rho_C = U1./a;
v_C = U2./U1;
T_C = (gamma-1)*(U3./U1 - (gamma/2).*v_C.^2);
P_C = rho_C.*T_C;
M_C = (v_C./sqrt(T_C));
Mass_flow_C = rho_C.*a.*v_C;
rho_throat(k) = rho_C(throat);
v_throat(k) = v_C(throat);
T_throat(k) = T_C(throat);
P_throat(k) = P_C(throat);
M_throat(k) = M_C(throat);
flow_rate_throat(k) = Mass_flow_C(throat);
figure(4);
if k==100
plot(x,Mass_flow_C,'m','linewidth',1)
hold on
else if k == 200
plot(x,Mass_flow_C,'y','linewidth',1)
hold on
elseif k == 500
plot(x,Mass_flow_C,'k','linewidth',1)
hold on
else if k == 1000
plot(x,Mass_flow_C,'w','linewidth',1)
hold on
else if k == 1200
plot(x,Mass_flow_C,'b','linewidth',1)
hold on
legend('100 Delta','200 Delta','500 Delta','1000 Delta','1200 Delta')
title('Non Dimensional mass flow rate distribution accross the nozzel at different time step - conservation form')
xlabel('Non Dimensional distance x/L')
ylabel('Mass flow rate in conservaitive form')
grid on
end
end
end
end
end
Time_C = toc
figure(5)
subplot(5,1,1)
plot(x,M_C,'m','linewidth',1)
ylabel({'mach number'})
xlabel({'macch number at throat'})
subplot(5,1,2)
plot(x,v_C,'y','linewidth',1)
ylabel({'Velocity'})
xlabel({'Velocity at throat'})
subplot(5,1,3)
plot(x,rho_C,'k','linewidth',1)
ylabel({'Density'})
xlabel({'Density at throat'})
subplot(5,1,4)
plot(x,T_C,'w','linewidth',1)
ylabel({'Temperature'})
xlabel({'Temperature at throat'})
subplot(5,1,5)
plot(x,P_C,'b','linewidth',1)
ylabel({'pressure'})
xlabel({'pressure at throat'})
figure(6)
plot(M_throat,'m','linewidth',1.5)
hold on
plot(rho_throat,'y','linewidth',1.5)
hold on
plot(T_throat,'k','linewidth',1.5)
hold on
plot(v_throat,'w','linewidth',1.5)
hold on
plot(P_throat,'b','linewidth',1.5)
legend('Mach Number','rho/rho_o','(v/a_o)','(T/T_o)','(P/P_o)')
xlabel('No.of time steps')
ylabel('M,rho/rho_o,v/a_o,T/T_o,P/P_o')
title('Timewise variation of the feild variables at throat')
grid on
endExplanation:
The purpose of the function is to calculate the values of
Mass flow rate, Pressure, Mach number, Density, velocity, Temperature, Density at the throat, the velocity at the throat, Temperature at the throat,
Mass flow rate at the throat, Pressure at the throat, Mach number at the throat by the known equations for all the grid points.
To determine the initial profile of Density & Temperature, for loop is used inputting the known conditions within an else-if command
Initial values of area, throat, velocity & Pressure are calculated using formulas.
In order to simplify and avoid errors, the equation is split into Solution Vector(U), Flux Vector(F) and Source Term().
• Now Outer Time loop is computed for loop to calculate all the grid points for various formulas inside the loop.
• MacCarmack technique is used to solve the numerical solution of governing equations such as Continuity, Momentum & Energy
Equation.
• The technique implies two methods for solving namely the Predictor method followed by the Corrector method.
• Predictor & Corrector method are assigned within the for loop individually for calculating values.
• The calculated values are averaged in the next step and updating the solution with old values.
• Boundary Conditions are applied within the for loop and subsequently used for calculating the Mass Flow Rate, Mach
Number & Pressure and also calculating the values throat section.
• Finally plotting is done for Distance(x) & Mass Flow Rate for different Timesteps.
• Program is ended.
Output :
Observation :
From the above plots we can conclude that initial variations observed in Non Conservative form for a certain period of time makes it consume higher amount of time to acheive stable solution when compared to Conservative form.
Variation :
Here by observing the above plots we could see that both the forms are in same phase giving the exact same solution for all the values.
Mass Flow rate :
Above plots shows that Mass Flow Rate at different Time steps along the Distance oscillates initially, however the solution becomes stable after the 500th Timestep for both the forms attaining a constant value.
Comparison :
Observation:
Eventhough the final solution of both the forms are almost similar, the Non Conservative form oscillates much when compared to Conservative form. The latter also gives almost stable solution along the length. By this we can conclude that Conservative form attains stability in lesser number of iterations with earliest time.
Grid Points
n = 31
|
Terms |
Exact_Sol |
Non_Conservative_Sol |
Conservative_Sol |
|
density |
0.634 |
0.63868 |
0.65203 |
|
Temperature |
0.833 |
0.83646 |
0.83893 |
|
pressure |
0.528 |
0.53423 |
0.54701 |
|
Mach number |
1 |
0.99938 |
0.97865 |
|
Velocity |
1 |
0.91401 |
0.89638 |
n= 41
|
Terms |
Exact_Sol |
Non_Conservative_Sol |
Conservative_Sol |
|
density |
0.634 |
0.63797 |
0.63362 |
|
Temperature |
0.833 |
0.8358 |
0.83019 |
|
pressure |
0.528 |
0.53322 |
0.52602 |
|
Mach number |
1 |
0.99964 |
1.0102 |
|
Velocity |
1 |
0.91389 |
0.92041 |
n=51
|
Terms |
Exact_Sol |
Non_Conservative_Sol |
Conservative_Sol |
|
density |
0.634 |
0.63765 |
0.63631 |
|
Temperature |
0.833 |
0.8355 |
0.832 |
|
pressure |
0.528 |
0.53276 |
0.52941 |
|
Mach number |
1 |
0.99981 |
1.0037 |
|
Velocity |
1 |
0.91388 |
0.91551 |
n = 91
|
Terms |
Exact_Sol |
Non_Conservative_Sol |
Conservative_Sol |
|
density |
0.634 |
0.6404 |
0.63793 |
|
Temperature |
0.833 |
0.83678 |
0.83332 |
|
pressure |
0.528 |
0.53587 |
0.53159 |
|
Mach number |
1 |
1.0028 |
1.0007 |
|
Velocity |
1 |
0.91729 |
0.91347 |
Conclusion :
Cumulating all the above results we can evidently conclude that Conservative form converges slightly faster than Non Conservative
form offering stable solution along the simulation. Yet both the forms are dependable in calculating the solutions which is almost precise.
From the above outputs, The following results are made
1)density,pressure and temperature at intial and throat values are known
2) 4 grid points are considered n=31,61,51,91 for grid dependence test
3) with increase in grid points solution is more accurate
4) Non conservative can be easily programmed because they are in the form (p.v,t)
5) in conservative form,equations are expressed in flux variables and modified to solution vectors
6) Non conservative is more accurate because primitive values are easily solved