Week 7 - Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
Any numerical solution of isentropic flow though a quasi ID subsonic-supersonic nozzle is unreasonable, while
analytic solution exists for these kind of problems. However we are trying here to establish the broadened capability
of CFD in various streams rather in representing the approximated form of anayltic solution. Here, we implement the
numerical solution with Conservative and Non-Conservative forms.
Governing Equations are marched in terms of PDE for both forms in order to use Time Marching solution for the
problem. These forms represent in numerical aspect only, if we analytic all the methods are same.
We consider the study, isentropic flow through a convergent divergent nozzle. The flow at the inlet to the
nozzle comes from reservoir. The cross section area of the reservoir is very large and hence velocity is very small. The
flow is locally subsonic the convergent section of the nozzle, Sonic at the throat and Supersonic at the divergent
section
Let us consider a small control volume of the nozzle as shown below
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 the nozzle, , 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 asuming uniform flow properties across any given
cross-section. Such flow is defined as Quasi one dimensional flow.
The governing continuty, momentum and energy equations in terms of non dimensional variables for this quassi one
dimensional staedy and isentropc flow can be expressed respectively as:
NON CONSEVATIVE FORM:-
CONSEVATIVE FORM:-
MacCarmack method:
To discretize the numerical solution we use the technique called MacCarmack method. The McCormack method is a variation of the two-step Lax-Wendroff Scheme but is much simpler in application. The application of McCormack method proceeds in two steps; a predictor step which is followed by a corrector step.
Let f(x) is equation to be solved. then the following steps are followed
i.Initial profile
Before the calculation, the initial profile for all the flow field variables which are velocity, density and temperature
has to be determined. These are guess values but should be specified with more awareness on the problem.
ii. Predictor step
(∂f/∂t)p=(f(i+1)-f(i))/Δx
iii. Solution update
fp=finitial+(∂f/∂t)*Δt
fp= flow field variables calculated after the predictor step
finitial = flow field variable in the initial step
Δt = time step size
iv. Corrector step
(∂fp/∂t)c = (f(i)-f(i-1))/Δx
(∂f/∂t)avg=((∂f/∂t)p+(∂fp/∂t)c)/2
v. Final solution update
fsol=finitial+(∂f/∂t)avg*Δt
Main program of solving Quasi 1-D subsonic-supersonic nozzle
%SOLVING QUASI ID SUPERSONIC NOZZLE FLOW EQUATIONS
%INPUT PARAMETERS
%No of grid points
n = 31;
gamma = 1.4;
x =linspace(0, 3, n) ;
dx = x(2)- x(1);
%No of time step
nt = 1400;
%courant Number
C = 0.5;
%FUNCTION FOR NON CONSERVATIVE FORM
tic;
[Mass_ flow_rate_NC, Pressure_NC, Mach_number NC, rho_NC, V_NC, T_NC, rho_throat _NC, v_throat_NC,T_throat _NC,Time_elapsed_NC = toc;
fprintf('Time Elasped for Non Conservative form: %0.4f', Time_eläpsed_NC);
hold on
%PLOTTING
%Timestep variation of properties at throat area in Non Conservative form
figure(2)
%Density
subplot (4,1,1)
plot(linspace(1,nt,nt), rho_throat_NC, 'color', 'm')
ylabel('Density')
legend ({'Density at throat'});
axis([0 1400 0 1.3]);
grid minor;
title('Timestep variation at throat area - NON CONSERVATIVE FORM')
%Temperature
subplot (4,1,2)
plot(linspace(1,nt,nt),T_throat_NC,'color','c')
ylabel('Temperature')
legend ({'Temperature at throat'});
axis([0 1400 0.6 1]);
grid minor;
%Pressure
subplot(4,1,3)
plot(linspace(1,nt,nt), Pressure_throat_NC, 'color', 'g')
ylabel('Pressure')
legend({'Pressure at throat'});
axis([0 1400 0.4 1.1]);
grid minor;
%Mach Number
subplot (4,1,4)
plot(linspace(1,nt,nt), Mach_number_throat _NC, 'color', 'y')
ylabel('Timesteps')
label('Mach Number')
legend({'Mach Number at throat'});
axis([0 1400 0.6 1.4]);
grid minor;
%steady state simulation for primitve values for Non conservative form
figure (3);
%Density
subplot(4,1,1);
plot(x,rho_NC,'color', 'm')
ylabe1('Density')
legend ({'Density'});
axis([0 3 0 1])
grid minor;
title('Variation of Flow Rate Distribution - NON CONSERVATIVE')
%Temperature
subplot(4,1,2);
plot(x, T_NC, 'color','c')
ylabel('Temperature')
legend({'Temperature'}) ;
axis([0 3 0 1])
grid minor;
%Pressure
subplot 4,1,3);
plot(x, Pressure_NC, 'color', 'g')
ylabel('Pressure')
legend({'Pressure'});
axis ([0 3 0 1])
grid minor;
%Mach number
subplot(4,1,4) ;
plot(x,Mach_number_NC, 'color','y')
xlabel('Nozzle Length')
ylabel('Mach number')
legend({'Mach number'});
axis ([0 3 0 4])
grid minor;
%FUNCTION FOR CONSERVATIVE FORM
tic;
[Mass_flow_rate_C, Pressure_C, Mach number_C, rho_C, v_C, T_C, rho_throat_C, v_throat_C, T_throat_C, Mass_flow_rate_C, Time_elasped_C = toc;
fprintf ('Time Elasped for Conservative form: %0.4f', Time_elasped_C);
%PLOTTING
%Timestep variation of properties at throat area in Conservative form
figure (5)
%Density
subplot(4,1,1)
plot(rho_throat_C,'color', 'm')
ylabel('Density')
legend({'Density at throat'});
axis([0 1400 0 1]);
grid minor;
title('Timestep variation at throat area - CONSERVATIVE FORM')
%Temperature
subplot(4,1,2)
plot(T_throat_C, 'color', 'c')
ylabel('Temperature')
legend ({'Temperature at throat'});
axis ([0 1400 0.6 1]) ;
grid minor;
%pressure
subplot(4,1,3)
plot(Pressure_throat_C,'color','g')
ylabel('Pressure')
legend({'Pressure at throat'}) ;
axis([o 1400 0.2 0.8]);
grid minor;
%Mach Number
subplot (4,1,4)
plot (Mach_number_throat_C, 'color', 'y')
xlabel ('Timesteps )
ylabel ('Mach Number')
legend ({'Mach Number at throat'});
axis([0 1400 0.6 1.4]) ;
grid minor;
%Steady state simulation for primitve values for Non conservative form
figure(6);
%Density
subplot(4,1,1);
plot(x,rho_C, 'color','m')
ylabel('Density')
legend ({'Density'});
axis([0 3 0 1])
grid minor;
title( Variation of Flow Rate Distribution - CONSERVATIVE FORM')
*Temperature
subplot(4,1,2) ;
plot(x,T_C, 'color', 'c')
ylabel' Temperature')
legend ({'Temperature'});
axis([0 3 0 1])
grid minor;
%Pressure
subplot(4,1,3) ;
plot(x, Pressure_C,'color', 'g')
ylabel('Pressure')
legend ({'Pressure });
axis([0 3 0 1])
grid minor;
%Mach number
subplot(4,1,4);
plot(x,Mach_number_C, 'color','y')
xlabel 'Nozzle Length')
ylabel ('Mach number')
legend ({'Mach number'});
axis([0 3 0 4])
grid minor;
%COMPARISON PLOT of Normalized Mass Flow Rate of both forms
figure(7);
hold on
plot(x,Mass_ flow_rate_NC,'r','linewidth',1.25)
hold on;
plot(x,Mass_flow_rate_C, 'b','linewidth', 1.25)
hold on;
grid on
legend ('Non Conservative', 'Conservative');
xlabel ('Length of the Domain')
label ('Mass Flow Rate')
title ('Comparison of Normalized Mass Flow Rate of both forms')Explanation:
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 toc command is used for calculating the time elasped 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 throat as well as Convergent/Divergent section.
Non conservative form
%Function for Non-Conservative form
Function [Mass _flow_rate_NC, Pressure_NC, Mach_number_NC,rho_NC,v_NC,T_NC,rho_throat_NC, v_throat_NC, T_throat_NC, Time_elapsed_NC,Mass _flow_rate_throat_NC, Pressure_throat_NC, Mach_number_throat_NC]
Intial boundary conditions
rho_NC = 1 - 0.3146*x;
T_NC = 1 -0.2341*x;
v_NC = (0.1+1.09*x).*T_NC.^0.5;
%Area & Throat profiles
a = 1 + 2.2*(x - 1.5).^2;
throat = find(a==1);
%Outer Time Loop
for k = 1:nt
%Time-step control
dt = min(C.*dx./(sqrt(T_NC)+v_NC));
%Saving a copy of old Values
rho_old = rho_NC;
T_old = T_NC;
v_old = v_NC;
%Predictor Method
for j = 2:n-1
%Simplifying the Equation's spatial terms
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;
dlog_a_dx = (log(a(j + 1)) - log(a(j)))/dx;
%continuity Equation
drho_dt_P(j) = -rho_NC(j)*dv_dx - rho_NC(j)*v_NC(j)*dlog_a_dx-v_NC(j)*drho_dx;
%Momentum Equation
dv_dt_P(j)=-v_NC(j)*dv_dx-(1/gamma)*(dT_dx + (T_NC(j)*drho_dx/rho_NC(j)));
%Energy Equation
dT_dt_P(j) = -v_NC(j)dT_dx-(gamma - 1)*T_NC(j)*(dv_dx + v_NC(j)*dlog_a_dx);
%Solution Update
rho_NC(j)= rho_NC(j)+drho_dt_P(j)*dt;
v_NC(j) = v_NC(j)+dv_dt_P(j)*dt;
T_NC(j) = T_NC(j)+dT_dt_P(j)*dt;
end
%corrector Method
for j = 2:n-1
%Simplifying the Equation's spatial terms
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;
dlog _a_dx = (log(a(j))-log(a(j - 1)))/dx;
%Continuity Equation
drho_dt_C(j) = -rho_NC(j)*dv_dx - rho_NC(j)*v_NC(j)*dlog_a_dx - v_NC(j)*drho_dx;
%Momentum Equation
dv_dt_C(j) = -v_NC(j)*dv_dx - (1/gamma)*(dT_dx + (T_NC(j)*drho_dx/rho_NC(j)));
%Energy Equation
dT_dt_C(j) = -v_NC(j)*dT_dx - (gamma- 1)*T_NC(j)*(dv_dx + v_NC(j)*dlog_a_dx);
end
%Compute Average Time derivative
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);
%Final Solution Update
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
%Applying Boundary Conditions
%At Inlet
rho_NC(1) = 1;
v_NC(1)= 2*v_NC(2)-v_NC(3);
T_NC(1) = 1;
%At outlet
rho_NC(n) = 2*rho_NC(n-1) - rho_NC(n-2) ;
v_NC(n) = 2*v_NC(n-1) - v_NC(n-2) ;
T_NC(n) = 2*T_NC(n-1) - T_NC(n-2);
%Defining Mass flow rate, Pressure and Mach number
Mass_flow_rate_NC = rho_NC.*a.*v_NC;
Pressure_NC = rho_NC.*T_NC;
Mach_number_NC = (v_NC./sqrt(T_NC));
%Calculating Variables at throat section
rho_throat_NC(k) = rho_NC(throat);
T_throat_NC(k) = T_NC(throat);
v_throat_NC(k) = v_NC(throat);
Mass_flow_rate_throat_NC(k) = Mass_flow_rate_NC(throat);
Pressure_throat_NC(k)= Pressure_NC(throat);
Mach_number_throat_NC(k) = Mach_number_NC(throat);
%PLOTTING
%Comparison of non-dimensional Mass flow rates at different time steps
figure(1);
if k==50
plot(x,Mass_flow_rate_NC,'m','linewidth',1.25);
hold on
elseif k==100
plot(x,Mass_flow_rate_NC,'c','linewidth',1.25);
hold on
elseif k==200
plot(x,Mass_flow_rate_NC,'g','linewidth',1.25);
hold on
elseif k==400
plot(x,Mass_flow_rate_NC,'y','linewidth' 1.25);
hold on
elseif k==800
plot(x,Mass_flow_rate_NC,'k','linewidth',1.25);
hold on
end
%Labelling
title('Mass flow rates at different Timesteps - NON CONSERVATIVE FORM')
xlabel('Distance')
ylabel('Mass flow rate')
legend ({'50^t^h Timestep', '100^t^h Timestep', '200^t^h Timestep', '400^t^h Timestep', '800^t^h Timestep'})
axis ([0 3 0 2])
grid on
end
endExplanation:
The purpose of the function is to calculate the values of Mass flow rate, Pressure, Mach number, Density, velocity, Temperature, Density at throat, velocity at throat,
Temperature at throat, Mass flow rate at throat, Pressure at throat, Mach number at 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 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 Predictor method followed by 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.
Conservative form
%Function for Conservative form
Function [Mass _flow_rate_NC, Pressure_NC, Mach_number_NC,rho_NC,v_NC,T_NC,rho_throat_NC, v_throat_NC, T_throat_NC, Time_elapsed_NC,Mass _flow_rate_throat_NC, Pressure_throat_NC, Mach_number_throat_NC]
%Determining Initial profile of Density and Temperature
for i = 1:n
if (x(i)>=0 && 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)
rho_C(i) = 0.634 - 0.3879*(x(i) - 1.5) ;
T_C(i) = 0.833 - 0.3507*(x(i) - 1.5) ;
end
end
%Area & Throat profiles
a = 1 + 2.2*(x - 1.5).^2;
throat = find(a==1);
%Initial profile for Velocity and Pressure
v_C = 0.59./(rho_C.*a) ;
Pressure_C = rho_C.*T_C;
%Calculate the Initial Solution Vectors(U)
U1 = rho_C.*a;
U2 = rho_C.*a.*v_C;
U3 = rho_C.*a.*(T_C./(gamma - 1)) + ((gamma/2).*v_C.^2);
%outer Time Loop
for k =1:nt
%Time-step control
dt = min(C.*dx./(sqrt(T_C)+v_C));
%Saving a copy of old Values
U1_old = Ul;
U2_old = U2:
U3_old = U3;
%Calculating the Flux vector(F)
F1 = U2;
F2 = ((U2.^2)./U1) + ((gamma - 1)/gamma)*(U3 - (gamma/2)*((U2.^2)./U1));
F3 = ((gamma*U2.*U3)./U1) - (gamma*(gamma - 1)*((U2.^3)./(2*U1.^2))) ;
%Predictor Method
for j = 2:n-1
%Calculating the Source term(J)
J2(j) = (1/gamma)*rho_C(j)*T_C(j)*((a(j + 1) - a(j))/dx);
%Continuity Equation
dU1_dt_P(j) = -((F1(j + 1) - F1(j))/dx);
%Momentum Equation
dU2_dt_P(j) = -((F2(j + 1) - F2(j))/dx) + J2(j);
%Energy Equation
dU3_dt_P(j) = -((F3(j + 1) - F3(j))/dx);
%updating New values
U1(j) = U1(j + dU1_dt_P(j)*dt;
U2(j) = U2(j)+ dU2_dt_P(j)*dt;
U3(j) = U3(j)+ dU3_dt_P(j)*dt;
end
%Updating Flux vectors
F1 = U2;
F2 = ((U2.^2)./U1) + ((gamma - 1)/gamma)*(U3 - ((gamma*0.5* (U2.^2))./(U1)));
F3 = ((gamma*U2.*U3)./U1) - ((gamma*0.5*(gamma - 1)*(U2.^3))./(U1.^2));
%Corrector Method
for j = 2:n-1
*updating the Source term(J)
J2(j)= (1/gamma)*rho_C(j)*T_C(j)*((a(j) - a(j - 1))/dx);
%simplifying the Equation's spatial terms
dU1-dt_C(j) = -((F1(j) - F1(1 - 1))/dx);
dU2_dt_C(j) = -((F2(j) - F2(j - 1))/dx) + J2(j) ;
dU3_dt_C(j) = -((F3(j) - F3(j - 1))/dx) ;
end
%Compute Average Time derivative
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) ;
%Final Solution Update
for i = 2:n-1
U1(i) = U1_old(i) + dU1_dt(i)*dt;
U2(i) = U2_old(i) + dU2_dt(i)*dt;
U3(i) = U3_old(i) + dU3_dt(i)*dt;
end
%Applying Boundary Conditions
%At Inlet
U1(1) = rho_C(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));
%At 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);
*Final Primitive value update
rho_C = U1./a;
v_C = U2./U1;
T_C = (gamma - 1)*((U3./U1 - (gamma/2)*(v_C).^2)) ;
%Defining Mass flow rate, Pressure and Mach number
Mass_flow_rate_C = rho_C.*a.*v_C;
Pressure_C = rho_C.*T_C;
Mach_number_C = (v_C./sqrt(T_C));
%Calculating Variables at throat section
rho_throat_C(k) = rho_C(throat);
T_throat_C(k) = T_C(throat);
v_throat_C(k) = v_C(throat);
Mass_flow_rate_throat_C(k) = Mass_flow_rate_C(throat);
Pressure_throat_C(k)= Pressure_C(throat);
Mach_number_throat_C(k) = Mach_number_C(throat);
%PLOTTING
%Comparison of non-dimensional Mass flow rates at different time steps
figure(4);
if k==50
plot(x,Mass_flow_rate_C,'m','linewidth',1.25);
hold on
elseif k==100
plot(x,Mass_flow_rate_C,'c','linewidth',1.25);
hold on
elseif k==200
plot(x,Mass_flow_rate_C,'g','linewidth',1.25);
hold on
elseif k==400
plot(x,Mass_flow_rate_C,'y','linewidth' 1.25);
hold on
elseif k==800
plot(x,Mass_flow_rate_C,'k','linewidth',1.25);
hold on
end
%Labelling
title('Mass flow rates at different Timesteps - CONSERVATIVE FORM')
xlabel('Distance')
ylabel('Mass flow rate')
legend ({'50^t^h Timestep', '100^t^h Timestep', '200^t^h Timestep', '400^t^h Timestep', '800^t^h Timestep'})
axis ([0 3 0.4 1])
grid on
end
endExplanation:
The purpose of the function is to calculate the values of Mass flow rate, Pressure, Mach number, Density, velocity, Temperature, Density at throat, velocity at throat,
Temperature at throat, Mass flow rate at throat, Pressure at throat, Mach number at throat by the known equations for all the grid points.
To determine initial profile of density and temperature, for loop is used inputing the unknown conditions withinan else-if command.
Initial values of area, throat, velocity & Pressure are calculated using formulas.
In order to simplify and avoid errors, equation is splitted into Solution Vector(U), Flux Vector(F) and Source Term (J).
Now Outer Time loop is computed for 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 Predictor method followed by 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.
Outputs:
Time steps variation of primitive variable at throat area
Non conservative form
Conservative form
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 of flow rate distribution:
Non conservative form
Conservative form
Observation:
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 at different timesteps:
Non conservative form
Conservative form
Observation:
Above plots shows that Mass Flow Rate at different Time steps along the Distance oscillates initially, however the solution becomes stable after the 400th Timestep for both the forms attaining a constant value.
Comparison plot of normalized mass flow rate of both Non censervative and conservative forms :
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.