Week 7 - Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
1. AIM
To Program and find the numerial solution of 1D supersonic nozzle folw simulation using Maccormack method. The following are the problem objectives :
- Find solution using Conservative and non-conservative forms
- 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
- Which method is faster?
2. INTRODUCTION
Let us consider a steady and isentropic flow throught the convergent diverent nozzle.
Inlet end has the reservoir where pressure and temperature are stagnant.
Velocity = 0 at inlet as CS of nozzle is large. Then the flow expands isendtropically and at the exit it reaches supersonic speed.
M<1 at convergent section; M=1 at throat of nozzle ; M>1 at divegent section
Governing equation for non concervative form:
- Continuity Equation
`(∂ ρ ')/ (∂ t ') = ρ( ' ∂ V )/(' ∂ x ' )− ρ ' V ' (∂ ( ln A ))/( ∂ x ' )− V '( ∂ p ')/( ∂ 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 ' )]`
Governing equation for Concervative form:
- Continuity Equation
`(∂ ( ρ ' A ' )) /(∂ t ' )+ (∂ ( ρ ' A ' V ' ) )/(∂ x ') = 0`
- Momentum Equation
`(∂ ( ρ ' A ' V ' ) )/(∂ t ') + (∂ ( ρ ' A ' V ' 2 + 1 /γ p ' A ' ) )/(∂ x ') = 1 /γ p ' (∂ A ')/( ∂ x ')`
- Energy Equation
`( ∂[ ρ ' ( (T ')/( γ − 1 )+ γ /2 V ' ^2 ) A ' ]) /(∂ t ') + (∂ [ ρ ' ( T ')/( γ − 1 )+ γ /2 V ' ^2 ) V ' A ' + p ' A ' V ' ] /(∂ x ' )= 0`
solving the Governing Equations
2.2 MATLAB CODE
Main Code - Solving 1D supersonic Nozzle
% MatLab Progra to simulate 1 supersonic Nozzle flow
clear all;
close all;
clc;
%Intitial Conditions
L = 3; %nozzle length
n = 31; % As in Anderson Text book
x = linspace(0,L,n);
dx = L/(n-1); %grid space
nt = 1400; % Time steps
throat = (n-1)/2; % Throat
gamma = 1.4; % Cp/Cv
A = 1 + 2.2*(x - 1.5).^2; % Area
C = 0.5; % Courant Number
%% Solution in conservative form
tic
[mc, Mc, Vc, Tc, Pc, rhoc, T_tc, rho_tc, P_tc, M_tc] = conservative_form(x,gamma,dx,A,n,nt,C,throat);
consv_form_time_taken = toc;
fprintf('\nComputation time by Conservative form: %0.4f',consv_form_time_taken);
%% Solution in non- Conservative form
tic
[m, M, V, T, P, rho, T_t, rho_t, P_t, M_t]= non_conservative_form(x,gamma,dx,A,n,nt,C,throat);
non_consv_form_time_taken = toc;
fprintf('\nComputation time by Non-Conservative form: %0.4f',non_consv_form_time_taken);
%% Plots
% convergence of Parameters at Nozzle throat at different time steps in Conservative and Non-Conservative Form:
figure(1);
subplot(4,1,1);
plot(1:nt,rho_tc,'b');
hold on
plot(1:nt,rho_t,'r');
title('Variations of density w.r.t Time');
legend('Conservative form', 'Non-Conservative form');
ylabel('rho/rho_{o}');
xlabel('Time steps');
axis([0 1400 0 1.3]);
subplot(4,1,2);
plot(1:nt,T_tc,'b');
hold on
plot(1:nt,T_t,'r');
title('Variations of Temperature w.r.t Time');
legend('Conservative form', 'Non-Conservative form');
ylabel('T/T_{o}');
xlabel('Time steps');
axis([0 1400 0.6 1]);
subplot(4,1,3);
plot(1:nt,P_tc,'b');
hold on
plot(P_t,'r');
title('Variations of Pressure w.r.t Time');
legend('Conservative form', 'Non-Conservative form');
ylabel('p/p_{o}');
xlabel('Time steps');
axis([0 1400 0.4 1.1]);
subplot(4,1,4);
plot(M_tc,'b');
hold on
plot(M_t,'r');
ylabel('M');
xlabel('Time steps');
title('Variations of Mach Number w.r.t Time');
legend('Conservative form','Non-Conservative form');
axis([0 1400 0.6 1.8]);
%% Variation of non dimensional flow parameters along x axis in Conservative and Non-Conservative Form:
figure(2);
subplot(4,1,1);
plot(x,Mc,'b');
hold on
plot(x,M,'*r');
title('Variation of non dimensional Mach number along x axis ');
legend('Conservative form', 'Non-Conservative form');
ylabel('M');
xlabel('Length(x)');
grid minor;
subplot(4,1,2);
plot(x,Pc,'b');
hold on
plot(x,P,'*r');
ylabel('p/p_{o}');
xlabel('Length(x)');
title('Variation of non dimensional Pressure ratio along x axis ');
legend('Conservative form', 'Non-Conservative form');
grid minor;
axis([0 3 0 1]);
subplot(4,1,3);
plot(x,rhoc,'b');
hold on
plot(x,rho,'*r');
ylabel('rho/rho_{o}');
title('Variation of non dimensional density ratio along x axis ');
legend('Conservative form', 'Non-Conservative form');
xlabel('Length(x)');
grid minor;
axis([0 3 0 1]);
subplot(4,1,4);
plot(x,Tc,'b');
hold on
plot(x,T,'*r');
ylabel('T/T_{o}');
xlabel('Length(x)');
title('Variation of non dimensional Temperature ratio along x axis ');
legend('Conservative form', 'Non-Conservative form');
grid minor;
axis([0 3 0 1]);
%% Comparison of mass flow rate quantity in Conservative and Non-Conservative form against Exact Form:
figure(3);
hold on;
plot(x, mc, 'r');
plot(x, m, 'b');
line([0 3], [0.579 0.579],'LineWidth', 1.5);
legend('Conservative form', 'Non-Conservative form', 'Exact form');
xlabel('x/L');
ylabel('rhoVA/rho_{o}V_{o}A_{0}');
title('non dimensional mass flow quantity in Conservative and Non-Conservative form against Exact form');
Sub code- Conservative Function
function [mc, Mc, Vc, Tc, Pc, rhoc, T_tc, rho_tc, P_tc, M_tc] = conservative_form(x,gamma,dx,A,n,nt,C,throat)
% Density and Temperature Initialize
for i = 1:n
if x(i) >= 0 && x(i) < 0.5
rhoc(i) = 1;
Tc(i) = 1;
elseif (x(i) >= 0.5 && x(i) < 1.5)
rhoc(i) = 1.0 -0.366*(x(i) - 0.5);
Tc(i) = 1.0 - 0.167*(x(i) - 0.5);
elseif (x(i) >= 1.5 && x(i) <= 3.0)
rhoc(i) = 0.634 - 0.3879*(x(i) - 1.5);
Tc(i) = 0.833 - 0.3507*(x(i) - 1.5);
end
end
% Velocity and Pressure Initialise
Vc = 0.59./(rhoc.*A);
Pc = rhoc.*Tc;
% solution vectors Parameters:
U1 = rhoc.*A;
U2 = rhoc.*A.*Vc;
U3 = rhoc.*A.*((Tc./(gamma - 1)) + (gamma/2)*Vc.^2);
% Time Iteration Loop
for k = 1:nt
U1_old = U1; % saving old solution vector
U2_old = U2;
U3_old = U3;
dt = min(C.*dx./(sqrt(Tc) + Vc));% optimum vaue of time step using courant number
% 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);
%% Predictor calculation:
for j = 2:n-1
% source term find J2:
J2(j) = (1/gamma).*rhoc(j).*Tc(j).*((A(j+1)-A(j))/dx);
%continuity equation
dU1dt_p(j) = -(F1(j+1) - F1(j))/dx;
%momentum equation
dU2dt_p(j) = -((F2(j+1) - F2(j))/dx) + J2(j);
%energy equation
dU3dt_p(j) = -(F3(j+1) - F3(j))/dx;
% Updating solution vectors:
U1(j) = U1(j) + dU1dt_p(j)*dt;
U2(j) = U2(j) + dU2dt_p(j)*dt;
U3(j) = U3(j) + dU3dt_p(j)*dt;
end
% Updating values of rho,V,T after predictor Step:
rhoc = U1./A;
Vc = U2./U1;
Tc = (gamma -1)*(U3./U1 - (gamma/2)*Vc.^2);
Pc = rhoc.*Tc;
%Updating flux - corrector step
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 method
for j = 2:n-1
% Find source term J2:
J2(j) = (1/gamma)*rhoc(j).*Tc(j)*((A(j)-A(j-1))/dx);
% Corrector step:
%continuity equation
dU1dt_c(j) = -(F1(j) - F1(j-1))/dx;
%momentum equation
dU2dt_c(j) = -((F2(j) - F2(j-1))/dx) + J2(j);
%energy equation
dU3dt_c(j) = -(F3(j) - F3(j-1))/dx;
end
% Finding the average of predictor and corrector steps:
dU1dt = 0.5*(dU1dt_p + dU1dt_c);
dU2dt = 0.5*(dU2dt_p + dU2dt_c);
dU3dt = 0.5*(dU3dt_p + dU3dt_c);
% Final solution vectors:
for j = 2:n-1
U1(j) = U1_old(j) + dU1dt(j)*dt;
U2(j) = U2_old(j) + dU2dt(j)*dt;
U3(j) = U3_old(j) + dU3dt(j)*dt;
end
% Applying the boundary conditions:
%Inlet
U1(1) = rhoc(1)*A(1);
U2(1) = 2*U2(2) - U2(3);
Vc(1) = U2(1)/U1(1);
U3(1) = U1(1)*((Tc(1)/(gamma-1))+((gamma/2)*Vc(1)^2));
% right side/outlet conditions:
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);
% primitive variables update:
rhoc = U1./A;
Vc = U2./U1;
Tc = (gamma -1)*(U3./U1 - (gamma/2)*Vc.^2);
Pc = rhoc.*Tc; % calculating Pressure
Mc = Vc./sqrt(Tc);% calculating Mach
mc = rhoc.*A.*Vc;% calculating mass
% Calculating values at the throat:
M_tc(k) = Mc(throat); % Mach Number
P_tc(k) = Pc(throat); % Pressure Ratio
T_tc(k) = Tc(throat); % Non-dimensional Temperature
rho_tc(k) = rhoc(throat); % Density Ratio
% Plotting mass flow rate at different time step:
figure(4);
hold on;
if k == 50
plot(x,mc,'b','linewidth',1.5);
elseif k == 100
plot(x,mc,'g','linewidth',1.5);
elseif k == 150
plot(x,mc,'r','linewidth',1.5);
elseif k == 200
plot(x,mc,'c','linewidth',1.5);
elseif k == 400
plot(x,mc,'m','linewidth',1.5);
elseif k == 700
plot(x,mc,'y','linewidth',1.5);
elseif k == 1400
plot(x,mc,'k','linewidth',1.5);
end
end
title('Comparison of non-dimensional mass flow rates at different time steps (Conservative Form)');
xlabel('x/L (Non-dimensional distance)');
ylabel('rhoVA/rho_{o}V_{o}A_{o}');
legend('50 Time steps','100 Time steps', '150 Time steps','200 Time steps','400 Time steps', '700 Time steps','1400 Time steps');
axis([0 3 0.4 1.2]);
endSub code- Non - conservative Function
function [m, M, V, T, P, rho, T_t, rho_t, P_t, M_t]= non_conservative_form(x,gamma,dx,A,n,nt,C,throat)
% Initial Condition
rho = 1-0.3146*x;
T = 1-0.2314*x;
V = (0.1+1.09*x).*T.^0.5;
dt = min(C.*dx./(sqrt(T)+V)); % optimum vaue of time step using courant number
% Time Loop
for k = 1:nt
rho_old = rho; % copy old values
T_old = T;
V_old = V;
% Predictor Step
for j = 2:n-1
dvdx = (V(j+1)-V(j))/dx;
drhodx = (rho(j+1) - rho(j))/dx;
dTdx = (T(j+1)-T(j))/dx;
dlog_a_dx = (log(A(j+1))-log(A(j)))/dx;
%Continuity Equation
drhodt_p(j) = -rho(j)*dvdx - rho(j)*V(j)*dlog_a_dx -V(j)*drhodx;
% Momentum Equation
dvdt_p(j) = -V(j)*dvdx - (1/gamma)*(dTdx + (T(j)/rho(j))*drhodx);
% Energy Equation
dTdt_p(j) = -V(j)*dTdx - (gamma-1)*T(j)*(dvdx + V(j)*dlog_a_dx);
% solution update - density , velocity and temperature
rho(j) = rho(j) + drhodt_p(j)*dt;
V(j) = V(j) + dvdt_p(j)*dt;
T(j) = T(j) + dTdt_p(j)*dt;
end
%corrector step
for j=2:n-1
dvdx = (V(j)-V(j-1))/dx;
drhodx = (rho(j) - rho(j-1))/dx;
dTdx = (T(j) - T(j-1))/dx;
dlog_a_dx = (log(A(j))-log(A(j-1)))/dx;
drhodt_c(j) = -rho(j)*dvdx - rho(j)*V(j)*dlog_a_dx - V(j)*drhodx;
dvdt_c(j) = -V(j)*dvdx - (1/gamma)*(dTdx + (T(j)/rho(j))*drhodx);
dTdt_c(j) = -V(j)*dTdx - (gamma-1)*T(j)*(dvdx + V(j)*dlog_a_dx);
end
% Calculating average of time derivative values
drhodt = 0.5*(drhodt_p + drhodt_c);
dvdt = 0.5*(dvdt_p + dvdt_c);
dTdt = 0.5*(dTdt_p + dTdt_c);
% Final solution
for i = 2:n-1
rho(i) = rho_old(i) + drhodt(i)*dt;
V(i) = V_old(i) + dvdt(i)*dt;
T(i) = T_old(i) + dTdt(i)*dt;
end
% Applying BC
% Inlet BC
V(1) = 2*V(2) - V(3);
rho(1) = 1;
T(1) = 1;
% Outlet BC
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);
% Non dimensional pressue, mass flow rate and mach number
P = rho.*T;
m = rho.*A.*V;
M = V./sqrt(T);
% Calcualtion of throat values
rho_t(k) = rho(throat);
T_t(k) = T(throat);
P_t(k) = P(throat);
M_t(k) = M(throat);
figure(5)
hold on;
if k == 50
plot(x,m,'b','linewidth',1.75);
hold on
grid on
elseif k==100
plot(x,m,'g','linewidth',1.5);
elseif k==150
plot(x,m,'r','linewidth',1.5);
elseif k==200
plot(x,m,'c','linewidth',1.5);
elseif k==400
plot(x,m,'m','linewidth',1.5);
elseif k==700
plot(x,m,'y','linewidth',1.5);
elseif k==1400
plot(x,m,'k','linewidth',1.5);
end
end
title('Comparison of non-dimensional mass flow rates at different time steps (Non-Conservative Form)');
xlabel('x/L (Non-dimensional distance through the nozzle)');
ylabel('rhoVA/rho_{o}V_{o}A_{o}');
legend('50 Time steps','100 Time steps', '150 Time steps','200 Time steps','400 Time steps', '700 Time steps','1400 Time steps');
axis([0 3 0 2]);
end
2.3 RESULTS AND DISCUSSION
1. Comparion of Conservative and Non Conservative form paramters -Variation of Density, Pressure,Velocity, Mach Number and Tempertature w.r.t length
Inference: It is clear from the results the profile matches with the one in the reference book and conservative and non conservative gives same profile.
2. Comparion of Conservative and Non Conservative form paramters -Density, Pressure,Velocity, Mach Number and Tempertature convergence at different time steps
Inference: The results obtained depicts the convergence of parameters for time steps 1400. Oscillations are visible initially in the parameters (density, temperature, pressure and mach number) under study due to the effect of expansion and compression that occurs in nozzle and gradually decreases with increase in time steps and reaches a stable condition. Conservative form show les time wise variations of quantities or disturbances when compared to non conservative form. The minimum no. of cycles for convergence would be 600 ie., after 600 iteratioons the paramters reach stable condition.
3. Comparison of non dimensional mass flow rates at different time steps for conservative and non conservative form
Inference: It is envisaged from the results that conservative form exhibit less variations in mass flow rate when compared to non conservative form despite of different initial condition. The conservative form results are closer to the exact value 0.579 at steady state. At about 700 time steps the values are converged for both conservative and non conservative form and closer to exact solution.
4. Comparion of normalized Conservative and Non Conservative form Mass Flow rate
Inference: It is depicted from the results that mass flow rate of conservative form is preserved and closer to the exact solution value 0.579 at steady state. The non conservative form exhibits oscillations at the inlet and oultet. Thus conservative results are better and more acceptable than non conservative results.
5. Grid Independent Test - Which Method is Fastest?
| Grid Independence Test | ||||||
| Grid Points (n) | Solution Form | ρ/ρo | P/Po | T/To | M | Computation Time (s) |
| 31 | Conservative | 0.6504 | 0.5463 | 0.8400 | 0.9818 | 2.6997 |
| 31 | Non-Conservative | 0.6387 | 0.5342 | 0.8365 | 0.9994 | 1.0514 |
| 61 | Conservative | 0.6383 | 0.533 | 0.8351 | 0.9993 | 1.8983 |
| 61 | Non-Conservative | 0.6377 | 0.5327 | 0.8354 | 0.9998 | 1.154 |
| NIL | Exact Solution | 0.634 | 0.528 | 0.833 | 1 | |
Inference: It is clear from the above table that even though when we increse the grid points from 31 to 61 or in other words reduce the grid size by half the change in the solution is visible only in 3rd and 4th significant digits. Most of the parameters values are same upto three significant digits. This confirms that after 31 grid points even though you increase the grid points or reduce grid size the solution is independent of grid size.
Seconldy it is evident that the Computation time of non conservative form is less compared to conservative form. So Non conservative form yield results faster than conservative form.
Minuimum No. of time steps simulation should run for convergence (for n=31 grid points) considering all quantities into account:
Conservative - 500 (approx)
Non conservative - 680 (approx)
700 time steps would be a conservative whole number value for convergence at steady state.
3. CONCLUSION
The program to solve the 1D supersonic Nozzle using conservative and non conservative form was developed and executed successfully.