CFD SOLUTION OF QUASI-ONE DIMENSIONAL SUBSONIC SUPERSONIC ISENTROPIC NOZZLE FLOW
Accomplishment in the project:
Numercial solutions are obtained for quasi-one dimensional subsonic-supersonic isentropic nozzle flow using a cfd technique named MacCormack's technique.
Heuristically speaking:
a). We are taking a flow problem with a known analytic solution which is a steady ,isentropic quasi one dimensional flow .
b). MacCormack's technique is employed and the final results that are obtained through this specific CFD technique are compared with the experimental results to determine the accuracy of the CFD solution.
c). And also to have a physical understanding of the flow field calculated using a CFD technique.
Quasi-one dimensional flow:
In a nozzle , the area changes as a function of x, hence in reality the flow field is minimum two dimensional.But by making the assumption that flow properties vary only with x, makes it equivalent to assuming uniform flow properties across any given section.This type of flow is defined as quasi-one dimensional flow.
Physics involved in the problem:
a). The flow is locally subsonic in the convergent section, sonic at the throat and supersonic at the divergent section.
b). The pressure ,temperature are considered to be as stagnation at the inlet,hence velocity is very small.And as flow expands isentropically to supersonic speeds machnumber increases(M>1).
c). For any mechanical system to accelerate particular amount of mass it needs force .Similarily in a nozzle to create a subsonic supersonic isentropic flow ,the force required comes with a specific pressure ratio across the nozzle.
Citation: J.J.Anderson, CFD The Basics with Applications,Chapter -7
1)Governing equations:(First eschelons of equations)
1.1. Conservation form of flow field governing equations for quasi one dimensional flow :(Non -dimensionalised form)
Continuity equation:
`(del(rho'*a'))/(delt')+(del(rho'a'v'))/(delx') = 0`
Momentum equation:
`(del(rho'*a'*v'))/(delt')+(del(rho'*a'*v'^2+(1/gamma)*p'*a'))/(delx')=1/gamma*p'(dela')/(delx')`
Energy equation:
`(del[rho'((e')/(gamma-1)+gamma/2*v'^2)*a'])/(delt')+(del[rho'((e')/(gamma-1)+gamma/2*v'^2)v'*a'+p'a'v'])/(delx')=0`
1.2.Non - Conservation form of flow field governing equations for quasi one dimensional flow:(Non dimensionalised form)
Continuity equation:
`(delrho')/(delt') = -rho'*(delv')/(delx')-rho'*v'*(del(lna'))/(delx')-v'(delrho')/(delx')`
Momentum equation:
`(delv')/(delt') = -v'*(delv')/(delx')-1/gamma*((delT')/(delx')+(T')/(rho')*(delrho')/(delx'))`
Energy equation:
`(delT')/(delt') = -v'*(delT')/(delx')-(gamma-1)*T'[(delv')/(delx')+v'*(del(lna'))/(delx')]`
Generic form for the above equations is given as:
`(delU_1)/(delt') = -(delF_1)/(delx')`
`(delU_2)/(delt') = -(delF_2)/(delx')+j_2`
`(delU_3)/(delt') = -(delF_3)/(delx')`
Where U is the solutions vector(contains the primitive varaiables), F is the flux vector and J_2 is the source term
The second eschelons of equations are the finite difference expressions that are implemented pertaining to maccormack's technique and put in the solvers code.
The third eschelons of equations are the boundary conditions and the time step calculation.
2). Matlab Syntax for solving and simulating the conservation and non conservation forms of the governing equations of quasi one dimensional subsonic supersonic isentropic flow through nozzle using the MacCormack method:
2.1. Main program:(for grid points n = 31 & n =61);
clear all
close all
clc
%%setting a command to check the grid dependency
grid_points = input('thirtyone,sixtyone')
%-------n=31--------%
if (grid_points == 31)
%%inputs
n = 31
x = linspace(0,3,n);
dx = x(2)-x(1);
gamma = 1.4;
%Area of the nozzle(non-dimensionalised)
a = 1+2.2*(x- 1.5).^2;
%analytical value of mass flow rate through nozzle
m_con_analytical = 0.579*ones(1,n);
%finding the throat
throat_index = find(a == 1);
%-----------------%%conservation form of governing equations-----------------%
tic
%%function representing the conservation form solver of flow field governing equations
[a k_con rho_con v_con t_con p_con m_con Mach_con rho_con_throat v_con_throat t_con_throat p_con_throat Mach_con_throat] = conservation_form_fluid_flow(n,x,dx,gamma,a,throat_index)
time_con = toc;
%%flow field values at the throat for different number of time steps
figure(2)
subplot(4,1,1)
plot(rho_con_throat,'r')
grid on
xlim([0 k_con])
ylim([0.5 1])
xlabel('number of iterations')
ylabel('rho/rho_{o}')
legend('density variation at throat')
title('Time wise variation of density at the throat(conservation form)')
subplot(4,1,2)
plot(t_con_throat,'b')
grid on
xlim([0 k_con])
ylim([0.6 1.0])
xlabel('number of iterations')
ylabel('t/t_{o}')
legend('temperature variation at throat')
title('Time wise variation of temperature at the throat(conservation form)')
subplot(4,1,3)
plot(p_con_throat,'k')
grid on
xlim([0 k_con])
ylim([0.4 1])
xlabel('number of iterations')
ylabel('p/p_{o}')
legend('pressure variation at throat')
title('Time wise variation of pressure at the throat(conservation form)')
subplot(4,1,4)
plot(Mach_con_throat,'g')
grid on
xlim([0 k_con])
ylim([0.9 1.3])
xlabel('number of iterations')
ylabel('M')
legend('Mach value variation at throat')
title('Time wise variation of Mach value at the throat(conservation form)')
%%steady state flow field values in conservation form across the nozzle
figure(3)
subplot(4,1,1)
plot(x,Mach_con,'g')
grid on
xticks([0 3])
yticks([0 1 Mach_con(end)])
xlabel('Distance through the nozzle')
ylabel('Mach number')
legend('Mach number in conservation form','location','Northwest')
title('Mach number across the nozzle')
subplot(4,1,2)
plot(x,t_con,'r')
grid on
xticks([0 3])
yticks([t_con(end) 0.833 1])
xlabel('distance through the nozzle')
ylabel('t/t_{o}')
legend('temperature in conservation form')
title('steady state temperature across the nozzle')
subplot(4,1,3)
plot(x,rho_con,'b')
grid on
xticks([0 3])
yticks([rho_con(end) 0.634 1])
xlabel('distance through the nozzle')
ylabel('rho/rho_{o}')
legend('density in conservation form')
title('steady state density across the nozzle')
subplot(4,1,4)
plot(x,p_con,'k')
grid on
xticks([0 3])
yticks([p_con(end) 0.528 1])
xlabel('distance through the nozzle')
ylabel('p/p_{o}')
legend('pressure in conservation form')
title('steady state pressure across the nozzle')
%-----------------%%non-conservation form of governing equations-----------------%
%%function representing the non conservation form solver of flow field governing equations
tic
[a k_nc rho v t p m_nc Mach rho_throat v_throat t_throat p_throat Mach_throat] = non_conservation_form_fluid_flow(n,x,dx,gamma,a,throat_index)
computational_time_non_con = toc;
%%flow field values at the throat for different number of time steps
figure(5)
subplot(4,1,1)
plot(rho_throat,'r')
grid on
xlim([0 k_con])
ylim([0.5 1])
xlabel('number of iterations')
ylabel('rho/rho_{o}')
legend('density variation at throat')
title('Time wise variation of density at the throat(non-conservation form)')
subplot(4,1,2)
plot(t_throat,'b')
grid on
xlim([0 k_con])
ylim([0.6 1.0])
xlabel('number of iterations')
ylabel('t/t_{o}')
legend('temperature variation at throat')
title('Time wise variation of temperature at the throat(non-conservation form)')
subplot(4,1,3)
plot(p_throat,'k')
grid on
xlim([0 k_con])
ylim([0.4 1])
xlabel('number of iterations')
ylabel('p/p_{o}')
legend('pressure variation at throat')
title('Time wise variation of pressure at the throat(non-conservation form)')
subplot(4,1,4)
plot(Mach_throat,'g')
grid on
xlim([0 k_con])
ylim([0.9 1.3])
xlabel('number of iterations')
ylabel('M')
legend('Mach value variation at throat')
title('Time wise variation of Mach at the throat(non-conservation form)')
%%steady state flow field values across the nozzle in non conservation form
figure(6)
subplot(4,1,1)
plot(x,Mach,'g')
grid on
xticks([0 3])
yticks([0 1 Mach(end)])
xlabel('Distance through the nozzle')
ylabel('Mach number')
legend('Mach number in non-conservation form','location','Northwest')
title('Mach number across the nozzle')
subplot(4,1,2)
plot(x,t,'r')
grid on
xticks([0 3])
yticks([t(end) 0.833 1])
xlabel('distance through the nozzle')
ylabel('t/t_{0}')
legend('temperature in non-conservation form')
title('steady state temperature across the nozzle')
subplot(4,1,3)
plot(x,rho,'b')
grid on
xticks([0 3])
yticks([rho(end) 0.634 1])
xlabel('distance through the nozzle')
ylabel('rho/rho_{0}')
legend('density in non-conservation form')
title('steady state density across the nozzle')
subplot(4,1,4)
plot(x,p,'k')
grid on
xticks([0 3])
yticks([p(end) 0.528 1])
xlabel('distance through the nozzle')
ylabel('p/p_{0}')
legend('pressure in non-conservation form')
title('steady state pressure across the nozzle')
%-------%%comaparison of mass flow rates in both conservation form and non conservation form------%
figure(7)
plot(x, m_con,x,m_nc,'r')
hold on
plot(x,m_con_analytical,'--')
grid on
xticks([0 0.3 0.6 0.9 1.2 1.5 1.8 2.1 2.4 2.7 3])
xlabel('distance through the nozzle')
ylabel('limits of steady mass flow rate')
legend('mass flow rate in conservation form','mass flow rate in non-conservation form','exact analytical solutiuon')
title('comparison of steady mass flow rate for conservation and non conservation forms across the nozzle')
fprintf('Minimum number of iterations in conservation form to reach steady state is %g',k_con)
fprintf('\nComputational time for conservation form solver is %g',time_con)
fprintf('\nMinimum number of iterations in non-conservation form to reach steady state is %g',k_nc)
fprintf('\nComputational time for non-conservation form solver is %g\n',computational_time_non_con)
end
%-------------n = 61----------------%
if (grid_points == 61)
%inputs
n = 61
x = linspace(0,3,n);
dx = x(2)-x(1);
gamma = 1.4;
a = 1+2.2*(x- 1.5).^2;
m_con_analytical = 0.579*ones(1,n);
throat_index = find(a==1);
%-----------------%%conservation form of governing equations-----------------%
tic
[a k_con rho_con v_con t_con p_con m_con Mach_con rho_con_throat v_con_throat t_con_throat p_con_throat Mach_con_throat] = conservation_form_fluid_flow(n,x,dx,gamma,a,throat_index)
time_con = toc;
%%Change in flow field values at the throat
figure(2)
subplot(4,1,1)
plot(rho_con_throat,'r')
grid on
xlim([0 k_con])
ylim([0.5 1])
xlabel('number of iterations')
ylabel('rho/rho_{o}')
legend('density variation at throat')
title('Time wise variation of density at the throat(conservation form)')
subplot(4,1,2)
plot(t_con_throat,'b')
grid on
xlim([0 k_con])
ylim([0.6 1.0])
xlabel('number of iterations')
ylabel('t/t_{o}')
legend('temperature variation at throat')
title('Time wise variation of temperature at the throat(conservation form)')
subplot(4,1,3)
plot(p_con_throat,'k')
grid on
xlim([0 k_con])
ylim([0.4 1])
xlabel('number of iterations')
ylabel('p/p_{o}')
legend('pressure variation at throat')
title('Time wise variation of pressure at the throat(conservation form)')
subplot(4,1,4)
plot(Mach_con_throat,'g')
grid on
xlim([0 k_con])
ylim([0.9 1.3])
xlabel('number of iterations')
ylabel('M')
legend('Mach variation at throat')
title('Time wise variation of Mach at the throat(conservation form)')
%%steady state flow properties across the nozzle in conservation form
figure(3)
subplot(4,1,1)
plot(x,Mach_con,'g')
grid on
xticks([0 3])
yticks([0 1 Mach_con(end)])
xlabel('Distance through the nozzle')
ylabel('Mach number')
legend('Mach number in conservation form','location','Northwest')
title('Mach number across the nozzle')
subplot(4,1,2)
plot(x,t_con,'r')
grid on
xticks([0 3])
yticks([t_con(end) 0.833 1])
xlabel('distance through the nozzle')
ylabel('t/t_{o}')
legend('temperature in conservation form')
title('steady state temperature across the nozzle')
subplot(4,1,3)
plot(x,rho_con,'b')
grid on
xticks([0 3])
yticks([rho_con(end) 0.634 1])
xlabel('distance through the nozzle')
ylabel('rho/rho_{o}')
legend('density in conservation form')
title('steady state density across the nozzle')
subplot(4,1,4)
plot(x,p_con,'k')
grid on
xticks([0 3])
yticks([p_con(end) 0.528 1])
xlabel('distance through the nozzle')
ylabel('p/p_{o}')
legend('pressure in conservation form')
title('steady state pressure across the nozzle')
%-----------------%%non-conservation form of governing equations-----------------%
tic
[a k_nc rho v t p m_nc Mach rho_throat v_throat t_throat p_throat Mach_throat] = non_conservation_form_fluid_flow(n,x,dx,gamma,a,throat_index)
computational_time_non_con = toc;
%%flow field values at the throat in non conservation form
figure(5)
subplot(4,1,1)
plot(rho_throat,'r')
grid on
xlim([0 k_con])
ylim([0.5 1])
xlabel('number of iterations')
ylabel('rho/rho_{o}')
legend('density variation at throat')
title('Time wise variation of density at the throat(non-conservation form)')
subplot(4,1,2)
plot(t_throat,'b')
grid on
xlim([0 k_con])
ylim([0.6 1.0])
xlabel('number of iterations')
ylabel('t/t_{o}')
legend('temperature non variation at throat')
title('Time wise variation of temperature at the throat(non-conservation form)')
subplot(4,1,3)
plot(p_throat,'k')
grid on
xlim([0 k_con])
ylim([0.4 1])
xlabel('number of iterations')
ylabel('p/p_{o}')
legend('pressure non variation at throat')
title('Time wise variation of pressure at the throat(non-conservation form)')
subplot(4,1,4)
plot(Mach_throat,'g')
grid on
xlim([0 k_con])
ylim([0.9 1.3])
xlabel('number of iterations')
ylabel('M')
legend('Mach variation at throat')
title('Time wise variation of Mach at the throat(non-conservation form)')
%%steady state flow field values across the nozzle in non-conservation form
figure(6)
subplot(4,1,1)
plot(x,Mach,'g')
grid on
xticks([0 3])
yticks([0 1 Mach(end)])
xlabel('Distance through the nozzle')
ylabel('Mach number')
legend('Mach number in non-conservation form','location','Northwest')
title('Mach number across the nozzle')
subplot(4,1,2)
plot(x,t,'r')
grid on
xticks([0 3])
yticks([t(end) 0.833 1])
xlabel('distance through the nozzle')
ylabel('t/t_{0}')
legend('temperature in non-conservation form')
title('steady state temperature across the nozzle')
subplot(4,1,3)
plot(x,rho,'b')
grid on
xticks([0 3])
yticks([rho(end) 0.634 1])
xlabel('distance through the nozzle')
ylabel('rho/rho_{0}')
legend('density in non-conservation form')
title('steady state density across the nozzle')
subplot(4,1,4)
plot(x,p,'k')
grid on
xticks([0 3])
yticks([p(end) 0.528 1])
xlabel('distance through the nozzle')
ylabel('p/p_{0}')
legend('pressure in non-conservation form')
title('steady state pressure across the nozzle')
%-------%%comaparison of mass flow rates in both conservation form and non conservation form------%
figure(7)
plot(x, m_con,x,m_nc,'r')
hold on
plot(x,m_con_analytical,'--')
grid on
xticks([0 0.3 0.6 0.9 1.2 1.5 1.8 2.1 2.4 2.7 3])
xlabel('distance through the nozzle')
ylabel('limits of steady mass flow rate')
legend('mass flow rate in conservation form','mass flow rate in non-conservation form','exact analytical solutiuon')
title('comparison of steady mass flow rate for conservation and non conservation forms across the nozzle')
fprintf('Minimum number of iterations in conservation form to reach steady state is %g',k_con)
fprintf('\nComputational time for conservation form solver is %g',time_con)
fprintf('\nMinimum number of iterations in non-conservation form to reach steady state is %g',k_nc)
fprintf('\nComputational time for non-conservation form solver is %g\n',computational_time_non_con)
end
2.2).Function to solve the conservation form of governing equations
function [a k_con rho_con v_con t_con p_con m_con Mach_con rho_con_throat v_con_throat t_con_throat p_con_throat Mach_con_throat] = conservation_form_fluid_flow(n,x,dx,gamma,a,throat_index)
%%variations in density ,temperature and velocity across the grid
for j = 1:n
if ((x(j)>=0) && (x(j)<=0.5))
rho_con(j) = 1;
t_con(j) = 1;
elseif((x(j)>=0.5) && (x(j)<=1.5))
rho_con(j) = 1 - 0.366*(x(j) - 0.5);
t_con(j) = 1 - 0.167*(x(j) - 0.5);
elseif ((x(j)>=1.5) && (x(j)<=3.0))
rho_con(j) = 0.634 - 0.3879*(x(j) - 1.5);
t_con(j) = 0.833 - 0.3507*(x(j) - 1.5);
end
end
v_con = 0.59./(rho_con.*a);
%%courant's number
C = 0.5;
%%synthesizing U_1,U_2 and U_3 by using variations in density,temperature and velocity
U_1 = rho_con.*a;
U_2 = rho_con.*a.*v_con;
U_3 = rho_con.*((t_con./(gamma-1)) + (gamma/2).*((v_con).^2)).*a;
%%convergence criterion
error_con = 1e9; %residual
tol = 1e-5; %tolerance
k_con = 0; %iteration counter
while (error_con>tol)
%iteration counter to establish the number of iterations it takes for the solver to reach staedy state flow
k_con= k_con+1 ;
m_con_old = rho_con.*a.*v_con;
dt = min(C*dx./(sqrt(t_con)+v_con));
U_1_OLD = U_1;
U_2_OLD = U_2;
U_3_OLD = U_3;
%%predictor step
%flux vectors defined in terms of dependant variables U_1,U_2,U_3
F_1 = U_2;
F_2 = ((U_2.^2)./(U_1)) + (((gamma-1)/gamma)*(U_3 - ((gamma/2)*((U_2.^2)./(U_1)))));
F_3 = ((gamma)*((U_2.*U_3)./(U_1))) - (((gamma*(gamma-1))/(2))*((U_2.^3)./(U_1.^2)));
for i = 2:n-1
%%source term
J_2(i) = ((1/gamma)*(rho_con(i))*(t_con(i)))*((a(i+1)-a(i))/dx);
%%continuity equation in generic form
du_1dt_p(i) = -(F_1(i+1)-F_1(i))/(dx);
%%momentum equation in generic form
du_2dt_p(i) = -((F_2(i+1)-F_2(i))/(dx))+J_2(i);
%%energy equationin its generic form
du_3dt_p(i) = -(F_3(i+1)-F_3(i))/(dx);
##Predicted values(barred quantities)
U_1(i) = U_1(i)+du_1dt_p(i)*dt;
U_2(i) = U_2(i)+du_2dt_p(i)*dt;
U_3(i) = U_3(i)+du_3dt_p(i)*dt;
end
%%decoding the predicted values of primitive variables
rho_con = U_1./a;
v_con = U_2./U_1;
t_con = (gamma-1)*((U_3./U_1)-((gamma/2)*(v_con).^2));
%%updating the barred quantities into the grid points
F_1 = U_2;
F_2 = ((U_2.^2)./(U_1)) + (((gamma-1)/gamma)*(U_3 - ((gamma/2)*((U_2.^2)./(U_1)))));
F_3 = ((gamma)*((U_2.*U_3)./(U_1))) - (((gamma*(gamma-1))/(2))*((U_2.^3)./(U_1.^2)));
%%corrector method
for c = 2:n-1
J_2(c) = ((1/gamma)*(rho_con(c))*(t_con(c)))*((a(c)-a(c-1))/(dx));
%%continuity equation in generic form
du_1dt_c(c) = -(F_1(c)-F_1(c-1))/(dx);
%%momentum equation in generic form
du_2dt_c(c) = -((F_2(c)-F_2(c-1))/(dx))+J_2(c);
%%energy equationin its generic form
du_3dt_c(c) = -(F_3(c)-F_3(c-1))/(dx);
end
%%Average derivatives
du_1dt_avg = 0.5*(du_1dt_p+du_1dt_c);
du_2dt_avg = 0.5*(du_2dt_p+du_2dt_c);
du_3dt_avg = 0.5*(du_3dt_p+du_3dt_c);
%%corrected values
for f = 2:n-1
U_1(f) = U_1_OLD(f)+du_1dt_avg(f)*dt;
U_2(f) = U_2_OLD(f)+du_2dt_avg(f)*dt;
U_3(f) = U_3_OLD(f)+du_3dt_avg(f)*dt;
end
%%Boundary conditions
%%Subsonic inflow boundary conditions
rho_con(1) = 1;
t_con(1) = 1;
U_1(1) = rho_con(1)*a(1);
U_2(1) = 2*U_2(2)-U_2(3);
v_con(1) = U_2(1)/U_1(1);
U_3(1) = U_1(1)*((t_con(1)/(gamma-1)) + (gamma/2)*v_con(1)^2);
%%supersonic outflow boundary conditions
U_1(end) = 2*(U_1(end-1))-(U_1(end-2));
U_2(end) = 2*(U_2(end-1))-(U_2(end-2));
U_3(end) = 2*(U_3(end-1))-(U_3(end-2));
%%Calculating the primitive variables(flow field values) from dependent variables(solutions vector elements)
rho_con = U_1./a;
v_con = U_2./U_1;
t_con = (gamma-1)*((U_3./U_1)-((gamma/2)*(v_con).^2));
%mass flow rate
m_con = rho_con.*a.*v_con;
%Mach number is given by the ratio of local velocity divided by the local speed of sound
Mach_con = v_con./sqrt(t_con);
%non-dimensionalized presssure
p_con = rho_con.*t_con;
%%convergence check or residual calculation with respect to time(number of iterations)
error_con = max(abs(m_con_old-m_con));
%%calculating the variation of flow field variable at the throat
rho_con_throat(k_con) = rho_con(throat_index);
v_con_throat(k_con) = v_con(throat_index);
t_con_throat(k_con) = t_con(throat_index);
p_con_throat(k_con) = p_con(throat_index);
Mach_con_throat(k_con) = Mach_con(throat_index);
%%variation in mass flow distribution across the nozzle at different time steps
H = figure(1);
hold on;
if k_con == 1
plot(x,m_con,'--r')
elseif k_con == 50
plot(x,m_con,'b')
elseif k_con == 100
plot(x,m_con,'y')
elseif k_con == 150
plot(x,m_con,'g')
elseif k_con == 200
plot(x,m_con,'m')
elseif k_con == 700
plot(x,m_con,'--k')
end
grid on
xticks([0 0.3 0.6 0.9 1.2 1.5 1.8 2.1 2.4 2.7 3.0])
yticks([0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8])
xlabel('non-dimensional distance through nozzle(x/L)')
ylabel('non-dimensional mass flow rate(rho*v*a)/(rho_{o}*v_{o}*a_{o})')
legend('1deltat','50deltat','100deltat','150deltat','200deltat','700deltat','location','northwest')
title('instantaneous distributions of non dimensional mass flow rate(conservation form)')
end
end 2.3). Function to solve the non conservation form of the governing equations:
function [a k_nc rho v t p m_nc Mach rho_throat v_throat t_throat p_throat Mach_throat] = non_conservation_form_fluid_flow(n,x,dx,gamma,a,throat_index)
%feeding the initial conditions
rho = 1 - 0.3146*x; %rho array
t = 1-0.2314*x; %t = temperature array
v = (0.1+1.09*x).*t.^(0.5); %velocity array
%courant's number
C = 0.5;
%%courant number based time step control for project
%%convergence criterion
error_nc = 1e9;
tol = 1e-5;
k_nc = 0;
%outer time loop
while (error_nc>tol)
k_nc =k_nc+1;
%%initial mass flow rate
m_nc_old = rho.*a.*v;
rho_old = rho;
v_old = v;
t_old =t;
%time step
dt = min(C*dx./(sqrt(t)+v));
%predictor method
for j =2:n-1
%continuity equation
% drho_dt_p = -rho(j)*(v(j+1)-v(j))/dx-rho(j)*v(j)*(log(a(j+1))-log(a(j)))/dx-v(j)*(rho(j+1)-rho(j)/dx)
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);
%Barred quantities
rho(j) = rho(j)+drhodt_p(j)*dt;
v(j) = v(j)+dvdt_p(j)*dt; %this way we are going to access the right values for velocity at the right nodes
t(j) = t(j)+dtdt_p(j)*dt;
end
%corrector method
for j =2:n-1
%%inserting hthe barred quantities into the grid points
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;
%continuity equation
drhodt_c(j) = -rho(j)*dvdx - rho(j)*v(j)*dlog_a_dx - v(j)*drhodx;
%momentum equation
dvdt_c(j) = -v(j)*dvdx - (1/gamma)*(dtdx + (t(j)/rho(j))*drhodx);
%energy equation
dtdt_c(j) = -v(j)*dtdx - (gamma-1)*t(j)*(dvdx+v(j)*dlog_a_dx);
end
%computing the average time derivative
drhodt_avg = 0.5*(drhodt_p+drhodt_c);
dvdt_avg = 0.5*(dvdt_p+dvdt_c);
dtdt_avg = 0.5*(dtdt_p+dtdt_c);
for i =2:n-1
%final solution update
rho(i) = rho_old(i)+drhodt_avg(i)*dt;
v(i) = v_old(i)+dvdt_avg(i)*dt;
t(i) = t_old(i)+dtdt_avg(i)*dt;
end
%apply the boundary conditions
%inlet
v(1) = 2*v(2)-v(3);
%outlet
v(end) = 2*v(end-1)-v(end-2);
rho(end) = 2*rho(end-1)-rho(end-2);
t(end) = 2*t(end-1)-t(end-2);
%Mach number is given by the ratio of local velocity divided by the local speed of sound
Mach = v./sqrt(t);
%%Non dimensional mass flow rate
m_nc = rho.*a.*v;
%non-dimensionalized presssure
p = rho.*t;
%%Convergence check or residual calculation with respect to time(number of iterations)
error_nc = max(abs(m_nc_old-m_nc));
%%calculating the variation of flow field variables at the throat
rho_throat(k_nc) = rho(throat_index);
v_throat(k_nc) = v(throat_index);
t_throat(k_nc) = t(throat_index);
p_throat(k_nc) = p(throat_index);
Mach_throat(k_nc) = Mach(throat_index);
L = figure(4);
hold on
%%variation in mass flow distribution across the nozzle at different time steps
if k_nc == 1
plot(x,m_nc,'--r')
elseif k_nc == 50
plot(x,m_nc,'b')
elseif k_nc == 100
plot(x,m_nc,'y')
elseif k_nc == 150
plot(x,m_nc,'g')
elseif k_nc == 200
plot(x,m_nc,'m')
elseif k_nc == 700
plot(x,m_nc,'--k')
end
grid on
xticks([0 0.3 0.6 0.9 1.2 1.5 1.8 2.1 2.4 2.7 3.0])
yticks([0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8])
xlabel('non-dimensional distance through nozzle(x/L)')
ylabel('non-dimensional mass flow rate(rho*v*a)/(rho_{o}*v_{o}*a_{o})')
legend('1deltat','50deltat','100deltat','150deltat','200deltat','700deltat')
title('instantaneous distributions of non dimensional mass flow rate(non-conservation form)')
end
end3)Final Numerical results :The steady state solution for conservation form:
3.1. Time wise variation of the density,temperature,pressure and mach number at the the nozzle throat:
Due to discretisation scheme that is employed for conservation form makes the plot appear smoother in the early times but number of time steps are higher to make it to the steady state solution.
3.2.) Instantaneous distributions of the non dimensional mass flow rate as a function of distance through the nozzle at six different times (approaching to steady state):
Variation of mass flow rate with increase in number of time steps is interesting to examine while it is approaching the steady state.The dotted line representing the initial conditions(0deltat) and as a procedure of time marching unsteady mass flow deviates from the initial conditions.
The line labeled 100(deltat) which is the result after 100 time steps shows a peak rise in mass flow rate across the nozzle but after 200 time steps it appears to be slowly a constant and after 700 iterations ,it is seems to be approaching the analytical value.
[when compared with the non conservation form plot the variations in mass flow distributions are pretty low which can be asserted is due to the assumed initial conditions]
3.3). Steady state flow field values:
4)Final Numerical results :The steady state solutions for non-conservation form:
4.1. Time wise variation of the density,temperature,pressure and mach number at the the nozzle throat:
The observation is that the at early time steps the changes are large which means the flow field values are driven by a stronger potential,hence changes are rapid.The exact values are almost collinear with the obtained steady state values.(From the citation)
4.2. Instantaneous mass flow distribution as a function of distance through the nozzle at different number of time steps:
Mass flow variations approach to the steady state attributed by the time marching procedure followed in solving the flow field values for different times.
The initial mass flow curve(dotted line) is quite distorted sine wave like,which is due to the product of the assumed initial values for density, velocity and parabolic variation of the nozzle area ratio.Radical changes can be seen from initial conditions(early times-0deltat) mass flow rate to the later times(100times) mass flow rate.
Mass flow variation keeps stumbling around inside the nozzle due to the unsteady nature of flow field values.And after 200 iterations it starts settling down to a constant value leading upto 700 time steps ,it appears a horizontal line(in red) across the graph indicating that convergence have been achieved.(density*velocity*area = constant)
4.3. The steady state flow field values across the nozzle
5)Comparison of normalized flow rate for bothe conservation and non conservation forms:
Here we can observe that the steady mass flow rate produced by conservative form is much more satisfactory due to its closeness to being a constant.where as on a magnified scale for both as done above in the plot,non-conservation form have oscillations at the beginning and end of the nozzle.
The steady state results are much closer to exact reults for conservation form when compared with non conservation form.Conservation form conserves mass throughout the flow field.
5).Grid Independence:(conservation form and non conservation form)
conditions at the nozzle throat
rho*/rho_0 t*/t_0 p*/p_0 M
Case 1:31 points 0.63859 0.83640 0.53411 0.99934
(Non-conservation form)
Case 2:61 points 0.63772 0.83547 0.53279 0.99983
(Non- conservation form)
Case 1:31 points 0.65062 0.84011 0.54659 0.98128
(Conservation form)
Case 2: 61 points 0.63851 0.83526 0.53332 0.99924
(Conservation form)
Experimental solution 0.634 0.833 0.528 1.000
The comparison between the steady state results for grid points 31 and 61 are made for both conservation and non-conservation forms between density,temperature,pressure and mach number,there is only a marginal improvement in the numerical solution when grid points are increased.
Hence we can summarise that the original calculations using grid points 31 is essentially grid independent for both the forms.The primitive variables are essentially closer to the exact values for non conservation form than conservation form.
6). Computational time and number of iterations to reach convergence:
The conservation form takes more number of timesteps(iterations) to reach steady state solution and second less computational time compared with non conservation form which is due to some time elapsed inside the non conservation solver for some unknown reasons.
(n=31)
Number of iterations(conservation form): 950
Computational time(conservation form): 5.97916 seconds
Number of iterations(Non -conservation form):767
Computational time(Non-conservation form): 6.58974 seconds
Conclusion:
1). Courant number is assumed to be 0.5 as it is giving closer values to the analytical solution from the author's experience.So other values(upto 1.1) are not effecting the numerical results and stability much till the number of time steps are adjusted for each of the courant value to make the same non dimensional time for different courant number.
2). Calculational effort is more in achieving solutions for conservation form rather than in achieving the solutions for non conservation forms.
3). There is no determination of the best amongst the both the slovers but non-conservation form results are distinctly closer to exact values.
4).Conservative form and non conservative form are equivalent in their derivative foms but when discretised they are not the same anymore.And more computational time in case of non conservative form can be attributed to the nature of the solver in which the number of terms are adding up everytime a ne w grid point is added and where as in conservative form they are cancelling out each other and being leftover only with boundary conditions.