Simulation of 1D Supersonic nozzle flow using MacCormack method   (MATLAB, Skill Lync, Hard)
Aim: To compute flow variables for 1D Subsonic Supersonic Flow Nozzle using MacCormack method for conservative and non-conservative equations.
Theory:
A subsonic-supersonic flow nozzle is a type of nozzle in which the inlet throat has Mach number of flow below 1, at throat it is equal to 1 and in divergent outlet region it is greater than 1. Here, it is assumed that flow parameters vary only in one direction(X-Axis). Although there exists a analytical solution for this 1D problem, here solution is obtained using MacCormack method. McCormack’s method is an explicit finite difference technique which is second order accurate in both space and time. For any problem, there are 2 types of equation forms, conservative and non-conservative. In this project both the forms are used and compared using output variables.
Main Code:
close all;
clear;
clc;
%% inputs
n = 31; %Number of nodes
nt = 1000; %Number of time steps
x = linspace(0,3,n);
%% Using conservative equations
[c_rho, c_temp, c_mfl, c_M, c_p, c_rho_th, c_v_th, c_temp_th, c_mfl_th, c_M_th, c_p_th, c_time] = conservative(n,nt,x);
%% Using Non-conservative equations
[nc_rho, nc_temp, nc_mfl, nc_M, nc_p, nc_rho_th, nc_v_th, nc_temp_th, nc_mfl_th, nc_M_th, nc_p_th, nc_time] = nonconservative(n,nt,x);
%% Data representaion
fprintf('The computation time for Conservative equations is: %f Secondsn',c_time);
fprintf('The computation time for Non-Conservative equations is: %f Secondsn',nc_time');
figh1 = figure(1);
figh1.WindowState = 'Maximized';
subplot(4,1,1)
plot(x,c_M,'r','LineWidth',1.5);
ylabel('M');
title('Flow Variables for Conservative equations');
grid on;
subplot(4,1,2)
plot(x,c_rho,'r','LineWidth',1.5);
ylabel('Density');
grid on;
subplot(4,1,3)
plot(x,c_temp,'r','LineWidth',1.5);
ylabel('Temperature');
grid on;
subplot(4,1,4)
plot(x,c_p,'r','LineWidth',1.5);
ylabel('Pressure');
xlabel('Distance from Inlet (x/L)');
grid on;
figh2 = figure(2);
figh2.WindowState = 'Maximized';
subplot(4,1,1)
plot(x,nc_M,'r','LineWidth',1.5);
ylabel('M');
title('Flow Variables for Non-conservative equations');
grid on;
subplot(4,1,2)
plot(x,nc_rho,'r','LineWidth',1.5);
ylabel('Density');
grid on;
subplot(4,1,3)
plot(x,nc_temp,'r','LineWidth',1.5);
ylabel('Temperature');
grid on;
subplot(4,1,4)
plot(x,nc_p,'r','LineWidth',1.5);
ylabel('Pressure');
xlabel('Distance from Inlet (x/L)');
grid on;
figh3 = figure(3);
figh3.WindowState = 'Maximized';
subplot(4,1,1)
plot(nc_M_th,'r','LineWidth',1.5);
hold on;
plot(c_M_th,'b','LineWidth',1.5);
ylabel('M');
legend('Non-conservative','Conservative');
title('Flow Variables for both equations at throat');
grid on;
subplot(4,1,2)
plot(nc_rho_th,'r','LineWidth',1.5);
hold on;
plot(c_rho_th,'b','LineWidth',1.5);
ylabel('Density');
grid on;
subplot(4,1,3)
plot(nc_temp_th,'r','LineWidth',1.5);
hold on;
plot(c_temp_th,'b','LineWidth',1.5);
ylabel('Temperature');
grid on;
subplot(4,1,4)
plot(nc_p_th,'r','LineWidth',1.5);
hold on;
plot(c_p_th,'b','LineWidth',1.5);
ylabel('Pressure');
xlabel('No. of iterations');
grid on;
figh4 = figure(4);
figh4.WindowState = 'Maximized';
for z = 1:nt
clf;
plot(x,nc_mfl(z,:),'r','LineWidth',1.5);
hold on;
plot(x,c_mfl(z,:),'b','LineWidth',1.5);
grid on;
legend('Non-conservative','Conservative');
title2 = sprintf('No. of iterations: %d',z);
title({'Mass Flow Rate in both types of equations';title2});
ylabel('Mass flow rate');
xlabel('Distance from Inlet (x/L)');
movieVector(z) = getframe(figh4);
end
%making video
myvid = VideoWriter('Mass Flow rate vs Time','MPEG-4');
myvid.FrameRate = 10;
open(myvid);
writeVideo(myvid, movieVector);
close(myvid);
Code & Theory for Conservative equations:
The non-dimensional governing equations are:
These equations are transformed by grouping variables:
The equations in the right are the condensed forms of equations. The 'U' variables are called solution vector variables, 'F' variables are called flux vector variables & 'J' is source vector variable. The variables like 'rho', 'e', 'v' are primitive variables. The values of 'F' & 'J' can be written in terms of 'U'.
Also primitive variables can be written as:
The initial conditions are:
Code:
function [rho,temp,mfl,M,p,rho_th,v_th,temp_th,mfl_th,M_th,p_th,time] = conservative(n,nt,x)
dx = x(2) - x(1);
thr_node = floor(((n-1))/2) + 1;
gamma = 1.4;
%% Initial profiles
a = 1 + (2.2*(x-1.5).^2); % area
for i = 1:n
if (x(i)>=0 && x(i)<=0.5)
rho(i) = 1;
temp(i) = 1;
elseif (x(i)>=0.5 && x(i)<=1.5)
rho(i) = 1 - 0.366*(x(i) - 0.5);
temp(i) = 1 - 0.167*(x(i) - 0.5);
elseif (x(i)>=1.5 && x(i)<=3.5)
rho(i) = 0.634 - 0.3879*(x(i) - 1.5);
temp(i) = 0.833 - 0.3507*(x(i) - 1.5);
end
end
v = 0.59./(rho.*a);
courant = 0.5; % courant number
dt = min(courant.*(dx./((temp.^0.5) + v))); % time step
%%MacCormack Method
%Defining initial solution vectors
u1 = rho.*a;
u2 = rho.*a.*v;
u3 = u1.*((temp./(gamma - 1)) + ((gamma/2).*(v.^2)));
tic
for k = 1:nt
u1_old = u1;
u2_old = u2;
u3_old = u3;
% flux terms defined in terms of u1,u2 and u3
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 Step
for j = 2:n-1
dadx = (a(j+1) - a(j))/dx;
j2 = rho(j)*temp(j)*dadx/gamma;
%Continuity Equation
du1dt_p(j) = -(f1(j+1) - f1(j))/dx;
%Momentum Equation
du2dt_p(j) = -(f2(j+1) - f2(j))/dx + j2;
%Energy Equation
du3dt_p(j) = -(f3(j+1) - f3(j))/dx;
%Solution Prediction
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
%Prediction of primitives
rho_p = u1./a;
v_p = u2./u1;
temp_p = (gamma - 1)*((u3./u1) - ((gamma/2)*(v_p.^2)));
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));
%Corrector Step
for j = 2:n-1
dadx = (a(j) - a(j-1))/dx;
j2 = rho_p(j)*temp_p(j)*dadx/gamma;
%Continuity Equation
du1dt_c(j) = -(f1(j) - f1(j-1))/dx;
%Momentum Equation
du2dt_c(j) = -(f2(j) - f2(j-1))/dx + j2;
%Energy Equation
du3dt_c(j) = -(f3(j) - f3(j-1))/dx;
end
%Average time derivative
du1dt = 0.5*(du1dt_p + du1dt_c);
du2dt = 0.5*(du2dt_p + du2dt_c);
du3dt = 0.5*(du3dt_p + du3dt_c);
%Final Solution
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
%Boundary Conditions
% inlet
u1(1) = rho_p(1)*a(1);
u2(1) = 2*u2(2) - u2(3);
v_p(1) = u2(1)/u1(1);
u3(1) = u1(1)*((temp_p(1)/(gamma-1)) + ((gamma/2)*v_p(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);
rho = u1./a;
v = u2./u1;
temp = (gamma - 1)*((u3./u1) - ((gamma/2)*(v.^2)));
mfl(k,:) = rho.*v.*a;
M = v./((temp).^0.5);
p = rho.*temp;
%Calculating values at throat
rho_th(k) = rho(thr_node);
v_th(k) = v(thr_node);
temp_th(k) = temp(thr_node);
mfl_th(k) = mfl(k,thr_node);
M_th(k) = M(thr_node);
p_th(k) = p(thr_node);
end
time = toc;
endCode & Theory for non-conservative equations:
The governing equations are:
These equations are non-dimensionalized:
For computation of these equations they are discritised by finite difference method:
These equations are solved by MacCormack Method. The time step used in both conservative and non-conservative approaches is calculated by finding values of time step at each node and then a minimum value of all time steps is choosen.
The Initial conditions are:
Code:
function [rho,temp,mfl,M,p,rho_th,v_th,temp_th,mfl_th,M_th,p_th,time] = nonconservative(n,nt,x)
dx = x(2) - x(1);
thr_node = floor(((n-1))/2) + 1;
gamma = 1.4;
%% Initial profiles
a = 1 + (2.2*(x-1.5).^2); % area
rho = 1-(0.3146*x); % density
temp = 1-(0.2314*x); % t = temperature
v = (0.1 + (1.09*x)).*temp.^0.5; % velocity
courant = 0.5; % courant number
dt = min(courant.*(dx./((temp.^0.5) + v))); % time step
%% MacCormack method
tic
for k = 1:nt
rho_old = rho;
v_old = v;
temp_old = temp;
%Predictor Step
for j = 2:(n-1)
dvdx = (v(j+1) - v(j))/dx;
drhodx = (rho(j+1) - rho(j))/dx;
dtdx = (temp(j+1) - temp(j))/dx;
dlogadx = (log(a(j+1)) - log(a(j)))/dx;
%Continuity equation
drhodt_p(j) = -(rho(j)*dvdx) - rho(j)*v(j)*dlogadx - v(j)*drhodx;
%Momentum equation
dvdt_p(j) = -v(j)*dvdx - (1/gamma)*(dtdx + (temp(j)/rho(j))*drhodx);
%Energy equation
dtdt_p(j) = -v(j)*dtdx - (gamma-1)*temp(j)*(dvdx + v(j)*dlogadx);
%predicted values
rho(j) = rho(j) + drhodt_p(j) * dt;
v(j) = v(j) + dvdt_p(j) * dt;
temp(j) = temp(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 = (temp(j) - temp(j-1))/dx;
dlogadx = (log(a(j)) - log(a(j-1)))/dx;
%Continuity equation
drhodt_c(j) = -(rho(j)*dvdx) - rho(j)*v(j)*dlogadx - v(j)*drhodx;
%Momentum equation
dvdt_c(j) = -v(j)*dvdx - (1/gamma)*(dtdx + (temp(j)/rho(j))*drhodx);
%Energy equation
dtdt_c(j) = -v(j)*dtdx - (gamma-1)*temp(j)*(dvdx + v(j)*dlogadx);
end
%Average Time Derivatives
drhodt = 0.5*(drhodt_p + drhodt_c);
dvdt = 0.5*(dvdt_p + dvdt_c);
dtdt = 0.5*(dtdt_p + dtdt_c);
%Corrected values
for j = 2:(n-1)
rho(j) = rho_old(j) + drhodt(j)*dt;
v(j) = v_old(j) + dvdt(j)*dt;
temp(j) = temp_old(j) + dtdt(j)*dt;
end
%Boundary conditions
%Inlet
v(1) = 2*v(2) - v(3);
%Outlet
rho(n) = 2*rho(n-1) - rho(n-2);
v(n) = 2*v(n-1) - v(n-2);
temp(n) = 2*temp(n-1) - temp(n-2);
p = rho.*temp;
M = v./(temp.^0.5);
mfl(k,:) = rho.*v.*a;
% Data required for representation purposes
drhodt_av(k) = abs(drhodt(thr_node));
dvdt_av(k) = abs(dvdt(thr_node));
rho_th(k) = rho(thr_node);
v_th(k) = v(thr_node);
temp_th(k) = temp(thr_node);
p_th(k) = p(thr_node);
M_th(k) = M(thr_node);
mfl_th(k) = mfl(k,thr_node);
end
time = toc;
endResults:
1)Variation of flow variables:
These are profiles of variables when they reach steady state. It can be seen that by using both forms of equations we get similar profiles of the variables. This is because we are solving the same problem using 2 different equations.
2) Variation of variables at throat with respect to time step
Here, we can compare the values obtained by both forms of equations:
In conservative form, the solution varies in small amounts with time step as compared to non-conservative form. The non-conservative equations give a wavy curve at beginning but as the time steps increase, it also becomes stable. The reason for this difference is that in conservation forms we conserve quantities. This conservation gives a steady state solution soon as compared to non-conservation.
3) Grid independence test and comparing the final results
This image shows that the program written is independent of grid points in mesh. Also it can be seen that non-conservative solutions give more accurate results than the conservative solution. This happens because the residuals that remain after each time step are reduced faster in non-conservative form. At last time step the residual was 10^(-6) in non-conservative form and was 10^(-3) in conservative form.
4) Command window output:
The computation time for Conservative equations is: 0.056445 Seconds
The computation time for Non-Conservative equations is: 0.029937 Seconds
Here, it can be seen that conservative equations take more time. This is because we have code and then decode the variables during each time step.
5 Variation of mass flow rate with time steps:
In this video we can see
1) Conservative solution does not vary in high amounts as compared to non-conservative solution.
2) Due to reason of residuals given above, it can be seen that conservative solution keeps on moving a little till the end.
3) The non-conservative solution at last time step seems to be shaped in U form but it is not wrong. If we compare it with the scale of Y axis, it can be seen that errors are small and profile can assumed to be straight for larger and practical scale.
Conclusion:
Both forms have their own merits and demerits and one should choose the form according to his requirements.