Nozzle flow using MacCormack s technique
The below link provides the Governing equations for both Conservation and Non-conservation forms of the nozzle flow using MacCormack\'s technique.
https://drive.google.com/open?id=1rpG3CAvFjnvN8XSDSL5cCHQVRUJbUe-g
Highlights of Report
1. Codes for both Conservation and Non-conservation form
2. Convergence study at nozzle throat (when area = 1)
3. Mass flow rate comparison for both the forms
4. Grid convergence study
1. Code snippet for Conservation and non-Conservation form
The code snippet is written for the grid spacing of 30 nodes and with the Courant no of 0.5. Time-marching method is used to obtain the solution. The final solution is represented after 1400 time steps.
a. Non-Conservation form
% Function for Predictor method
function [continuity,momentum,energy] = non_conservation_p(v,rho,t,a,j,dx,gamma)
dvdx = (v(j+1)-v(j))/dx;
drhodx = (rho(j+1)-rho(j))/dx;
dtdx = (t(j+1)-t(j))/dx;
dlog_adx = (log(a(j+1))-log(a(j)))/dx;
%Continuty equation
continuity = -rho(j)*dvdx-rho(j)*v(j)*dlog_adx - v(j)*drhodx;
%Momentum equation
momentum = -v(j)*dvdx - (1/gamma)*(dtdx + (t(j)/rho(j))*drhodx);
%Energy equation
energy = -v(j)*dtdx - (gamma-1)*t(j)*(dvdx + v(j)*dlog_adx);
end
% Function for Corrector method
function [continuity,momentum,energy] = non_conservation_c(v,rho,t,a,j,dx,gamma)
dvdx = (v(j)-v(j-1))/dx;
drhodx = (rho(j)-rho(j-1))/dx;
dtdx = (t(j)-t(j-1))/dx;
dlog_adx = (log(a(j))-log(a(j-1)))/dx;
%Continuty equation
continuity = -rho(j)*dvdx-rho(j)*v(j)*dlog_adx - v(j)*drhodx;
%Momentum equation
momentum = -v(j)*dvdx - (1/gamma)*(dtdx + (t(j)/rho(j))*drhodx);
%Energy equation
energy = -v(j)*dtdx - (gamma-1)*t(j)*(dvdx + v(j)*dlog_adx);
end
%% Main programme
%% Numerical Solution to Quasi-1D Non-Conservation form Nozzle flow
%% using MacCormack\'s technique
clear all
close all
clc
% Inputs
% n = input(\'No of grid points: \');
n = 31;
L = 3;
x = linspace(0,L,n);
dx = L/(n-1);
gamma = 1.4;
% Initial profiles
rho = 1-0.3146*x;
t = 1-0.2314*x; % t = temp
v = (0.1+1.09*x).*t.^0.5;
% Area
a = 1+2.2*(x - 1.5).^2;
% Time steps
nt = 1400;
% C = input(\'Input the required Courant no: \');
C = 0.5;
dt = min((C.*dx)./(sqrt(t)+v));
% Outer time loop
for k = 1:nt
rho_old = rho;
v_old = v;
t_old = t;
%% Predictor method
for j = 2:n-1
% Predictor method Governing equations
[drho_dt_p(j),dvdt_p(j),dtdt_p(j)] = non_conservation_p(v,rho,t,a,j,dx,gamma);
%Solution update
rho(j) = rho(j) + drho_dt_p(j)*dt;
v(j) = v(j) + dvdt_p(j)*dt;
t(j) = t(j) + dtdt_p(j)*dt;
end
%% Corrector method
for j = 2:n-1
% Corrector method Governing equations
[drho_dt_c(j),dvdt_c(j),dtdt_c(j)] = non_conservation_c(v,rho,t,a,j,dx,gamma);
end
% Compute the average time derivative
dvdt = 0.5*(dvdt_p + dvdt_c);
drhodt = 0.5*(drho_dt_p + drho_dt_c);
dtdt = 0.5*(dtdt_p+dtdt_c);
% Solution update
for i = 2:n-1
v(i) = v_old(i) + dvdt(i)*dt;
rho(i) = rho_old(i) + drhodt(i)*dt;
t(i) = t_old(i) + dtdt(i)*dt;
end
% Boundary conditions
% Inlet B.C.
v(1) = 2*v(2) - v(3);
%Outlet B.C.
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);
m = rho.*v.*a;
mach = v./(sqrt(t));
P = rho.*t;
% Solution at throat (area = 1)
rho_th(k) = rho(16);
temp_th(k) = t(16);
Pr_th(k) = P(16);
mach_th(k) = mach(16);
v_th(k) = v(16);
m_th(k) = m(16);
k_th = 1:nt;
end
% figure(1)
% plot(k_th,m_th,\'-m\'); hold on;
% legend(\'Non-conservation\');
figure
subplot(2,2,1);
plot(x,v); grid on;
xlabel(\'x/L\');
ylabel(\'v/a_o\');
title(\'Velocity profile\');
subplot(2,2,2);
plot(x,mach); grid on;
xlabel(\'x/L\');
ylabel(\'Mach no m\');
title(\'Mach no profile\');
subplot(2,2,3);
plot(x,rho); grid on;
xlabel(\'x/L\');
ylabel(\'{\\rho_e}/{\\rho_o}\');
title(\'Density profile\');
subplot(2,2,4);
plot(x,m); grid on;
xlabel(\'x/L\');
ylabel(\'{\\rho}av/{\\rho_o}a_ov_o\');
title(\'Mass flow rate profile\');
suptitle(\'Parameter profile for 1D nozzle flow for Non-Conservation form\');
%% Writing output in word file
fid = fopen(\'nonconservation.doc \',\'w\');
fprintf(fid,\'Parameter for 1D Conservation form nozzle flow \\n\');
fprintf(fid,\'for grid size %d and Courant no %f \\n\\n\', n,C);
fprintf(fid,\'Velocity Density Temp Mach no Massflowrate \\n\\n\');
fprintf(fid, \'%f %f %f %f %f \\n\', [v;rho;t;mach;m]);
fclose(fid);
The different parameters (Mass, density, Mach no, Velocity) variation for the non-conservation form is shown below. The parameters are non-dimensionalized and are varied with respect to the length of the nozzle.
Fig 1: Non-conservation flow parameter variation
b. Conservation form
% Function for Preditor value
function [continuity,momentum,energy] = conservation_p(t,a,i,dx,gamma,f1,f2,f3,rho)
dadx = (a(i+1)-a(i))/dx;
df1dx = (f1(i+1)-f1(i))/dx;
df2dx = (f2(i+1)-f2(i))/dx;
df3dx = (f3(i+1)-f3(i))/dx;
dlog_adx = (log(a(i+1))-log(a(i)))/dx;
j2 = (rho.*t.*dadx)./gamma;
% Continuity equation
continuity = -df1dx;
% Momentum equation
momentum = -df2dx + j2(i);
%Energy equation
energy = -df3dx;
end
% Function for Corrector value
function [continuity,momentum,energy] = conservation_c(t,a,i,dx,gamma,f1,f2,f3,rho)
dadx = (a(i)-a(i-1))/dx;
df1dx = (f1(i)-f1(i-1))/dx;
df2dx = (f2(i)-f2(i-1))/dx;
df3dx = (f3(i)-f3(i-1))/dx;
dlog_adx = (log(a(i))-log(a(i-1)))/dx;
j2 = (rho.*t.*dadx)./gamma;
% Continuity equation
continuity = -df1dx;
% Momentum equation
momentum = -df2dx + j2(i);
%Energy equation
energy = -df3dx;
end
%% Main Function
% Solving to Quasi 1D nozzle flow using MacCormak technique
% with the usage of Conservation form of Governing equation
clear all
close all
clc
% Initialization
% n = input(\'No of grid points: \');
n = 31;
L = 3;
x = linspace(0,L,n);
dx = L/(n-1);
gamma = 1.4;
%% Initial profiles and Conditions
% Area
a = 1+2.2*(x - 1.5).^2;
%Density and temp Initial conditions
rho = zeros(1,n);
t = zeros(1,n);
for i = 1:n
if (x(i)>=0 && x(i)<=0.5)
rho(i) = 1;
t(i) = 1;
elseif (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);
elseif (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
% Velocity condition
v = 0.59./(rho.*a);
% Initial conditions for Solution vectors
e = t;
p = rho.*t;
u1 = rho.*a;
u2 = rho.*a.*v;
u3 = rho.*a.*((e/(gamma-1))+((gamma/2)*v.^2));
% Time steps
nt = 1400;
%C = input(\'Input the required Courant no: \');
C = 0.5;
dt = min((C.*dx)./(sqrt(t)+v));
%% Outer time loop
for k = 1:nt
u1_old = u1;
u2_old = u2;
u3_old = u3;
% Flux terms in pure form
f1 = u2;
term1 = (u3./u1) - ((gamma/2)*((u2./u1).^2));
f2 = ((u2.^2)./u1)+(u1.*term1.*(gamma-1))./(gamma);
term2 = (gamma*(gamma-1)*(u2.^3))./(2*u1.^2);
f3 = ((gamma*u2.*u3)./u1)-term2;
%% Predictor Method
for i = 2:n-1
% Governing equations
[du1_dt_p(i),du2_dt_p(i),du3_dt_p(i)] = conservation_p(t,a,i,dx,gamma,f1,f2,f3,rho);
%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
rho = u1./a;
v = u2./u1;
t = (gamma-1)*((u3./u1)-((gamma*v.^2)/2));
% Update flux terms in pure form
f1 = u2;
term1 = (u3./u1) - ((gamma/2)*((u2./u1).^2));
f2 = ((u2.^2)./u1)+(u1.*term1.*(gamma-1))./(gamma);
term2 = (gamma*(gamma-1)*(u2.^3))./(2*u1.^2);
f3 = ((gamma*u2.*u3)./u1)-term2;
% Corrector Method
for i = 2:n-1
% Governing equations
[du1_dt_c(i),du2_dt_c(i),du3_dt_c(i)] = conservation_c(t,a,i,dx,gamma,f1,f2,f3,rho);
end
% Compute the average time derivative
du1dt = 0.5*(du1_dt_p + du1_dt_c);
du2dt = 0.5*(du2_dt_p + du2_dt_c);
du3dt = 0.5*(du3_dt_p + du3_dt_c);
% Solution update
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
% Apply the boundary conditions
% Inlet B.C.
u1(1) = rho(1)*a(1);
u2(1) = 2*u2(2) - u2(3);
v(1) = u2(1)/u1(1);
u3(1) = u1(1)*((t(1)/(gamma-1))+((gamma/2)*v(1)^2));
%Outlet B.C.
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);
% Obtaining primitives value i.e. rho, v, t and p
rho = u1./a;
v = u2./u1;
t = (gamma-1)*((u3./u1)-((gamma*v.^2)/2));
p = rho.*t;
mach = v./(sqrt(t));
mass = rho.*v.*a;
% Solution at throat (area = 1)
rho_th(k) = rho(16);
temp_th(k) = t(16);
Pr_th(k) = p(16);
mach_th(k) = mach(16);
v_th(k) = v(16);
mass_th(k) = mass(16);
k_th = 1:nt;
end
% plot(k_th,mass_th,\'-b\');
figure
subplot(2,2,1);
plot(x,v); grid on;
xlabel(\'x/L\');
ylabel(\'v/a_o\');
title(\'Velocity profile\');
subplot(2,2,2);
plot(x,mach); grid on;
xlabel(\'x/L\');
ylabel(\'Mach no m\');
title(\'Mach no profile\');
subplot(2,2,3);
plot(x,rho); grid on;
xlabel(\'x/L\');
ylabel(\'{\\rho_e}/{\\rho_o}\');
title(\'Density profile\');
subplot(2,2,4);
plot(x,mass); grid on;
xlabel(\'x/L\');
ylabel(\'{\\rho}av/{\\rho_o}a_ov_o\');
title(\'Mass flow rate profile\');
suptitle(\'Parameter profile for 1D nozzle flow for Non-Conservation form\');
%% Writing output in word file
fid = fopen(\'conservation.doc \',\'w\');
fprintf(fid,\'Parameter for 1D Conservation form nozzle flow \\n\');
fprintf(fid,\'for grid size %d and Courant no %f \\n\\n\', n,C);
fprintf(fid,\'Velocity Density Temp Mach no Massflowrate \\n\\n\');
fprintf(fid, \'%f %f %f %f %f \\n\', [v;rho;t;mach;mass]);
fclose(fid);
The different parameters (Mass, density, Mach no, Velocity) variation for the conservation form is shown below. The parameters are non-dimensionalized and are varied with respect to the length of the nozzle.
Fig 2: Conservation form flow parameter variation
2. Convergence study at the nozzle throat
The convergence study was on both the forms based on the no of time steps at the throat region where the area is one. The solution was represented in terms of plots which are shown below:
Fig 3: Convergence for different Parameters with time
Fig 3 represents the Convergence profile for different parameters with respect to time. The following observations can be made
a. At early times, the time derivatives are large and hence the values have oscillated. These oscillations are associated with various unsteady compression and expansion waves which propagate through the nozzle during the transient process.
b. At later times, the time derivatives rapidly grow small, changing by six orders of magnitude over a time span of 1400 time steps. Since we have used a finite no of time steps there will be some little amount of residual which is neglected.
c. From the mach no Convergence, we can clearly see that at initial time steps, the Mach no is not equal to 1 and after some time steps, it tends towards the value of 1 which represents the Converged solution.
d. From the above profiles, we can clearly say that the minimum no of simulations for the solution to get into convergence is around 700~800 time steps. However, in order to get more accuracy, 1000 time steps would be ideal.
3. Mass flow rate Comparison:
Fig 4: Mass flow rate Comparison
a. Fig 4 represents the mass flow rate comparison for both Conservation and Non-conservation forms.
b. From the fig, it can be observed that the conservation form gives a distribution that is much closer to being a constant. Whereas, the non-conservation form has a magnitude of higher range with some surplus oscillations at both the inflow and outflow boundaries.
c. It can be concluded that the Conservation form does a better job of preserving mass throughout the flow field, mainly because the mass flow itself is one of the dependent variables in the equations and the mass flow is a primary result from these equations.
d. With respect to Non-conservation form, the dependent variables in the non-conservation form of the equations are primitive variables, and the mass flow is obtained only as a secondary result.
Faster method
The computational effort is much higher in Conservation form than in non-conservation form. This is because it needs to decode the primitive variables from the flux variables. But such decoding is not necessary while solving non-conservation form. So, based on the complexity, the time required to solve the conservation form is more than the non-conservation form.
4. Grid convergence study
The grid convergence study in done for 3 different grids with the nodes involving 31,51 and 71 for the Courant no 0.5. The study is done for mass flow as the other parameters does not show much variation. The below fig represents the grid convergence for non-conservation form.
Fig 5: Grid dependency on the Mass flow rate for non-Conservation form
From the grid dependency study, we can observe that as the no of grids is increased, the steady-state numerical solution are slightly closer to the constant values. But care should be taken with respect to Courant no also. This is the case for both Conservation and Non-conservation form. This solution seems to get closer to exact solution as the grid are increased as per John D. Anderson textbook.
Conclusion
1. The conservation form yields a better mass flow distribution. Also, the conservation form does a better job of conserving the mass.
2. The non-conservation form is faster since the complexity of solving the conservation form is more.
3. As the no of grids is increased, with ideal courant no, the solution tends towards the exact solution.