Flow Through Convergent-Divergent Nozzle
Flow Through Nozzle
Introduction to Nozzle
A nozzle is a device designed to control the direction or characteristics of a fluid flow (to increase velocity) as it exists (or enters) an enclosed chamber or pipe. A nozzle is often a pipe or tube of varying cross sectional area, and it can be used to direct or modify the flow of a fluid(liquid or gas). Nozzles are frequently used to control the rate of flow, speed, direction, mass, shape, and pressure of the stream that emerges from them. In a nozzle, the velocity of the fluid increases at the expense of its pressure energy.
Flow through nozzles occurs in a variety of engineering applications:
- Convergent nozzles are used for subsonic and sonic flows. They can also be used as flow measuring and flow regulating devices.
- Convergent-Divergent nozzles are used for Supersonic flows. It is comparitively difficult to obtain convergent-divergent passages in compressors and turbine blade rows.
Here we are interested in Convergent-Divergent nozzle as we are writing a code to solve the 1D Supersonic nozzle flow equations using Macormack Method.
Figure 4.11 shows the flow from an infinite reservoir through a convergent-divergent nozzle to an exhaust chamber. Though the main subject here is isentropic flow,it will be observed that flow through a convergent-divergent nozzle is not isentropic under all condition.The static pressure distributions for various values of the pressure ratio `p_(b)/p_(0)` are shown in curves a,b,c,...,h,i and j.
In curves a and b the pressure ratio `p_e/p_0`,across the nozzle is such that the flow is accelerating only up to the throat;the diverging part acts as a diffuser through which the pressure rises to the back pressure `p_b`. In curve critical conditions(`p_t=p*`;`M_t`=1) are reached at the throat but the back pressure is such that the diverging part still acts as a diffuser.
Curve h corresponds to the design value of the back pressure.All other curves are for off-
design conditions.Here the flow is supersonic in the entire divergent part of the nozzle.The flow in the curves a,b,c,...,h,i and j is isentropic but non-isentropic in other curves.
After the conditions of curve c, when the back pressure is further lowered expansion takes
place to supersonic velocity beyond the throat to a point where a discontinuity in the flow occurs; this is on account of the\"off design”\"value of the back pressure and the supersonic flow as in curves d and e.In these cases at some point the minimum pressure reached in supersonic expansion has to rise to the back pressure. But it is known that for the
deceleration of supersonic flow the passage should be convergent. Therefore, the existing shape of the nozzle downstream of this point is incompatible with the required process;as a result of this the flow readjusts itself to the shape of the passage by suddenly becoming subsonic. The pressure,temperature and density suddenly rise to values which are compatible with the subsonic flow immediately downstream of the supersonic flow region.Such a sudden change of supersonic flow to subsonic occurs through a plane of discontinuity between the supersonic and subsonic flow regions.This plane of discontinuity is called a shock wave. The flow through the shock wave is no longer isentropic.
When the back pressure is lowered further the shock wave moves downstream till it reaches the exit as in curves f and g;here the back pressure is still higher than the design value.The nozzle exit pressure `p_e`,rises to the back pressure `p_b`,through a shock wave outside the nozzle exit section as shown in curves f and g.
Curves i and j correspond to back pressure lower than the nozzle exit pressure.The conditions in the nozzle remain unchanged and the adjustment of the nozzle exit pressure to the exhaust chamber pressure occurs through an expansion wave outside the nozzle.
Highlights
In this project we are using Macormack Method to solve the above case of Convergent Divergent Nozzle equation. In computational fluid dynamics, the MacCormack method is a widely used discretization scheme for the numerical solution of hyperbolic partial differential equations. This method essentially uses a predictor step and a corrector step which will be depicted in the code below. This method provides an accuracy of second order without actually using second order derivative and there by reducing the complexity of the process.
Here we are asked to study the fluid properties as it passes through the convergent-divergent nozzle. We are asked to track the properties in Conservative and Non-Conservative types of flows. Though both conservative and non-conservative are mathematically equal, they are different in numerical descritisation techniques. One form can be converted to another exactly.
Conservative Flow
Consider a fluid flow through a pathway or pipe. Consider a small control volume inside the pipe. In conservative flow we consider this small control volume to be stationery and we find the properties at this control volume. The conservative form of flow is more accurate than the other counterpart. It also works well with discontinuties(shock) in the domain.
The main governing equations are:
Continuity Equation : `(del.rhoA)/(delt) + (del(rhoVA))/(delx) = 0`
Momentum Equation :`(del(rhoAV))/(delt) + (del(rhoAV^2 + (1/gamma)P.A))/(delx) = (1/gamma)*P*(delA)/(delx)`
Energy Equation : `(del[rho(T/(gamma-1)+(gamma)/2.V^2)A])/(delt) + (del[rho(T/(gamma-1)+(gamma)/2*V^2)*V*A +P*A*V])/(delt)=0`
The above equations have been furthur simplified in the John D Anderson textbook page number 339. These steps eases the use of the above equations in programming. Let us define the elements of the solutions vector U, the flux vector F, and the source term J as follows:
The term e' stands for non dimensional internal energy that is `e' = e/e_0`, where `e_0 = c_vT_0 = (RT_0)/(gamma-1)`.
The term 'A' stands for area, 'V' for velocity, 'p' for pressure and `rho`for density.
Now with these above expressions we can rewrite the Continuity, Momentum and Energy Equations respectively.
In conservative form , the dependent variable (the variable for which we directly obtain numbers) are not primitive variable . To obtain primitive variables (`rho`,V,T,p) we must decode the elements `U_1,U_2 and U_3` as follows
Note that `e' = T'`as we have substituted in the above equations.
Flux terms
As per the John D Anderson textbook page number 341 & 342, we have the flux terms as:
`F_1 = U_1`
`F_2 = (U_2^2)/U_1 + (gamma-1)/gamma * (U_3-(gamma*U_2^2)/(2*U_1))`
`F_3 = gamma*(U_2*U_3)/U_1 - (gamma(gamma-1)*U_2^3)/(2*U_1^2`
`J_2 = (gamma-1)/gamma * (U_3 - (gamma*U_2^2)/(2*U_1))*(del(ln A'))/(delx'`
The flow properties at the downstream, supersonic outflow boundary are obtained by linear extrapolation from the two adjacent internal points. Let 'N' be the grid point at the outflow boundary then:
`(U_1)_N = 2(U_1)_(N-1)-(U_1)_(N-2)`
`(U_2)_N = 2(U_2)_(N-1)-(U_2)_(N-2)`
`(U_3)_N = 2(U_3)_(N-1)-(U_3)_(N-2)`
Initial conditions
Non-Conservative Flow
Consider a fluid flow through a pathway or pipe. Consider a small control volume inside the pipe. In non-conservative flow we consider this small control volume to be moving and compute the properties.
The main governing equations are:
Continuity Equation : `(delrho)/(delt) = - rho*(delV)/(delx)-rho*V(dellnA)/(dx) - V*(delrho)/(delx)`
Momentum Equation : `(delV)/(delt) = -V*(delV)/(delx) - (1/(gamma))*((delT)/(delx) + (T/P*(delrho)/(delx)))`
Energy Equation : `(delT)/(delt) = -V*(delT)/(delx) - (gamma-1)*T*((delV)/(delx)+V*(rholnA)/(delx))`
Initial conditions
Area = `1+2.2(x-1.5)^2, 0<=x<=3`
Assumptions
- The flow is considered to be Quasi-One dimensional flow as the properties aren\'t varied across y axis. Only variation along x axis is considered.
- A good amount of pressure ratio is maintained between chocked pressure and back pressure to allow expansion of the flow and avoid shock waves.
- The flow is assumed to be isentropic.
- Friction and heat transfer are not considered for simplicity.
Grid Dependency Test
For performing grid dependency test, first we set the value for n as 31 and run the program. We will obtain the respective graphs for n = 31. Later on we change the value of n to 61 and run the program again. The graphs obtained from the code where n = 61 is pasted below in the Outputs setion. These graphs are compared and suitable inference is reached.
MATLAB Main Code
%Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
clear all
close all
clc
%Inputs
n = 31; % change to 61 for grid dependency test
x = linspace(0,3,n); %Mesh Generation
dx = x(2)-x(1); %dist bw 2 points
gamma = 1.4;
nt = 1000; %time steps
%node of throat
throat_node = ((n-1)/2)+1;
%input your solver
%solver = 1;
%Conservative Function
%solver == 1;
tic
[rho_cons,v_cons,T_cons,mass_cons,p_cons,mach_cons,rho_throat1,p_throat1,T_throat1,mach_throat1,v_throat1,mass_throat1] = conservation_sec(x,dx,n,gamma,nt,throat_node);
time_cons = toc;
fprintf(\'computation time of conservative form = %d\',time_cons);
%Non Conservative Function
%solver == 2;
tic
[rho_non_cons,v_non_cons,T_non_cons,mass_non_cons,p_non_cons,mach_non_cons,rho_throat2,p_throat2,T_throat2,mach_throat2,v_throat2,mass_throat2] = nonconservation(x,dx,n,gamma,nt,throat_node);
time_non_cons = toc;
fprintf(\'computation time of Non Conservative form = %d\',time_non_cons);
%Plotting for Throat
figure(3)
p1 = subplot(4,1,1);
plot(linspace(1,nt,nt),rho_throat1);
hold on
plot(linspace(1,nt,nt),rho_throat2,\'color\',\'r\');
hold on
ylabel(\'Density\', \'fontsize\',14);
legend(p1,{\'Conservative form\' , \'Non-Conservative form\'});
%axis([0 1500]);
title(\'Time variation at throat\', \'fontsize\',14);
p2 = subplot(4,1,2)
plot(linspace(1,nt,nt),T_throat1);
hold on
plot(linspace(1,nt,nt),T_throat2,\'color\',\'r\');
hold on
ylabel(\'Temperature\',\'fontsize\',14);
legend(p2,{\'Conservative form\' , \'Non-Conservative form\'});
%axis([0 1500]);
p3 = subplot(4,1,3);
plot(linspace(1,nt,nt),p_throat1);
hold on
plot(linspace(1,nt,nt),p_throat2,\'color\',\'r\');
hold on
ylabel(\'Pressure\',\'fontsize\',14);
legend(p3,{\'Conservative form\' , \'Non-Conservative form\'});
%axis ([0 1500]);
p4 = subplot(4,1,4);
plot(linspace(1,nt,nt),mach_throat1);
hold on
plot(linspace(1,nt,nt),mach_throat2,\'color\',\'r\');
hold on
ylabel(\'Mach\',\'fontsize\',14);
legend(p4,{\'Conservative form\' , \'Non-Conservative form\'});
%axis([0 1500]);
%plotting for x axis
figure(4)
q1 = subplot(4,1,1);
plot(x,p_cons);
hold on
plot(x,p_non_cons,\'color\',\'r\',\'markersize\',4);
grid on
ylabel(\'Pressure\',\'fontsize\',14);
xlabel(\'Length of Nozzle\',\'fontsize\',8);
legend(q1,{\'Conservative form\' , \'Non-Conservative form\'});
title(\'Flow variation\', \'fontsize\',14);
q2 = subplot(4,1,2);
plot(x,mach_cons);
hold on
plot(x,mach_non_cons,\'color\',\'r\',\'markersize\',4);
grid on
ylabel(\'Mach\',\'fontsize\',14);
xlabel(\'Length of Nozzle\',\'fontsize\',8);
legend(q2,{\'Conservative form\' , \'Non-Conservative form\'});
q3 = subplot(4,1,3);
plot(x,T_cons);
hold on
plot(x,T_non_cons,\'color\',\'r\',\'markersize\',4);
grid on
ylabel(\'Temperature\',\'fontsize\',14);
xlabel(\'Length of Nozzle\',\'fontsize\',8);
legend(q3,{\'Conservative form\' , \'Non-Conservative form\'});
q4 = subplot(4,1,4);
plot(x,rho_cons);
hold on
plot(x,rho_non_cons,\'color\',\'r\',\'markersize\',4);
grid on
ylabel(\'Density\',\'fontsize\',14);
xlabel(\'Length of Nozzle\',\'fontsize\',8);
legend(q4,{\'Conservative form\' , \'Non-Conservative form\'});
%Mass flow rate comparison
figure(5)
plot(x,mass_cons);
hold on
plot(x,mass_non_cons,\'color\',\'r\',\'markersize\',4);
hold on
ylabel(\'Mass flow rate\',\'fontsize\',14);
xlabel(\'Distance\',\'fontsize\',14);
legend({\'Conservative form\' , \'Non-Conservative form\'});
title(\'Mass Flow Rate Comparison\', \'fontsize\',14);
grid on
Conservative Function Code
function [rho_cons,v_cons,T_cons,mass_cons,p_cons,mach_cons,rho_throat1,p_throat1,T_throat1,mach_throat1,v_throat1,mass_throat1] = conservation_sec(x,dx,n,gamma,nt,throat_node)
%CFL number
CFL = 0.5;
%area subjected
ar = 1+2.2*(x-1.5).^2;
%initial conditions
for j = 1:n;
if(x(j)>=0 && x(j)<=0.5)
rho_cons(j) = 1;
T_cons(j) = 1;
elseif (x(j)>=0.5 && x(j)<=1.5)
rho_cons(j) = 1-0.366.*(x(j)-0.5);
T_cons(j) = 1-0.167.*(x(j)-0.5);
elseif (x(j)>=1.5 && x(j)<=3.5)
rho_cons(j) = 0.634-0.3879.*(x(j)-1.5);
T_cons(j) = 0.833-0.3507.*(x(j)-1.5);
end
end
v_cons = 0.59./(rho_cons.*ar);
%Time step
%dt = min(CFL.*(dx./((T_cons.^0.5) + v_cons)));
%solution vector
U1 = rho_cons.*ar;
U2 = rho_cons.*ar.*v_cons;
U3 = U1.*((T_cons./(gamma-1)) + ((gamma/2).*(v_cons.^2)));
%Time loop
for k = 1:nt
%Time step
dt = min(CFL.*(dx./((T_cons.^0.5) + v_cons)));
%Updating the old values of U
U1_old = U1;
U2_old = U2;
U3_old = U3;
%Flux vector as in page 339 of Anderson text book
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
%replace the spatial derivatives with forward differencing
for h = 2:n-1
df1_dx(h) = (F1(h+1)-F1(h))/dx;
df2_dx(h) = (F2(h+1)-F2(h))/dx;
df3_dx(h) = (F3(h+1)-F3(h))/dx;
dlog_a_dx(h) = (log(ar(h+1))-log(ar(h)))/dx;
dar_dx(h) = (ar(h+1)-ar(h))/dx;
%Source vector
J2(h) = (1/gamma)*(rho_cons(h)*T_cons(h))* dar_dx(h);
%continuity equation
dU1_dt_p(h) = -df1_dx(h);
%momentum equation
dU2_dt_p(h) = -df2_dx(h)+J2(h);
%energy equation
dU3_dt_p(h) = -df3_dx(h);
%updating values
U1(h) = U1(h) + (dU1_dt_p(h))*dt;
U2(h) = U2(h) + (dU2_dt_p(h))*dt;
U3(h) = U3(h) + (dU3_dt_p(h))*dt;
end
rho_cons = U1./ar;
v_cons = U2./U1;
T_cons = (gamma-1)*((U3./U1)-((gamma*v_cons.^2)/2));
F1 = U2;
F2 = ((U2.^2)./U1) + (((gamma-1)/gamma)*(U3-(gamma*U2.^2)./(2*U1)));
F3 = ((gamma*U2.*U3)./U1)-(gamma*(gamma-1)*(U2.^3))./(2*U1.^2);
%corrector step
for p = 2:n-1
df1_dx(p) = (F1(p)-F1(p-1))/dx;
df2_dx(p) = (F2(p)-F2(p-1))/dx;
df3_dx(p) = (F3(p)-F3(p-1))/dx;
dlog_a_dx(p) = (log(ar(p))-log(ar(p-1)))/dx;
dar_dx(p) = (ar(p)-ar(p-1))/dx;
%Source vector
J3(p) = (1/gamma)*(rho_cons(p)*T_cons(p))*dar_dx(p);
%continuity equation
dU1_dt_c(p) = -df1_dx(p);
%momentum equation
dU2_dt_c(p) = -df2_dx(p)+J3(p);
%energy equation
dU3_dt_c(p) = -df3_dx(p);
end
%average values
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) = dU1_dt(i)*dt + U1_old(i);
U2(i) = dU2_dt(i)*dt + U2_old(i);
U3(i) = dU3_dt(i)*dt + U3_old(i);
end
%inlet boundary conditions
U1(1) = rho_cons(1)*ar(1);
U2(1) = 2*U2(2)-U2(3);
v_cons(1) = U2(1)/U1(1);
U3(1) = U1(1)*((T_cons(1)/(gamma-1))+((gamma/2)*v_cons(1)^2));
%outlet boundary 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);
%updateing old values
rho_cons = U1./ar;
v_cons = U2./U1;
T_cons = (gamma-1)*((U3./U1) - ((gamma*v_cons.^2)/2));
%Mach number
mach_cons = v_cons./(T_cons.^(0.5));
%Mass flow rate
mass_cons = rho_cons.*ar.*v_cons;
%Pressure
p_cons = rho_cons.*T_cons;
%Throat conditions
%rho_throat1,p_throat1,T_throat1,mach_throat1,v_throat1,mass_throat1
rho_throat1(k) = rho_cons(throat_node);
p_throat1(k) = p_cons(throat_node);
T_throat1(k) = T_cons(throat_node);
mach_throat1(k) = mach_cons(throat_node);
v_throat1(k) = v_cons(throat_node);
mass_throat1(k) = mass_cons(throat_node);
%To compare mass flow rate at different time steps
figure(1)
hold on
if k == 50
plot(x,mass_cons,\'b\',\'linewidth\',2)
elseif k == 100
plot(x,mass_cons,\'r\',\'linewidth\',2)
elseif k == 150
plot(x,mass_cons,\'c\',\'linewidth\',2)
elseif k == 200
plot(x,mass_cons,\'g\',\'linewidth\',2)
elseif k == 700
plot(x,mass_cons,\'y\',\'linewidth\',2)
elseif k == 1000
plot(x,mass_cons,\'.\',\'color\',\'k\',\'markersize\',10)
title(\'Mass flow rate conservative\');
xlabel(\'Non dimensional distance through nozzle\')
ylabel(\'Non dimensional mass flow\')
legend(\'50dt\',\'100dt\',\'150dt\',\'200dt\',\'700dt\',\'1000dt\')
end
end
end
Non-Conservative Function Code
function [rho_non_cons,v_non_cons,T_non_cons,mass_non_cons,p_non_cons,mach_non_cons,rho_throat2,p_throat2,T_throat2,mach_throat2,v_throat2,mass_throat2] = nonconservation(x,dx,n,gamma,nt,throat_node);
%CFL number
CFL = 0.5;
%initial condition
rho_non_cons = 1 - 0.3146 * x; %Density
T_non_cons = 1 - 0.2314 * x; %Temperature
v_non_cons = (0.1+(1.09*x)).*(sqrt(T_non_cons)); %Velocity
%Area
ar = (1 + 2.2 * ((x-1.5).^2));
%TIME LOOP
for k = 1:nt
%dt = min(abs(CFL.*(dx./((T_non_cons.^0.5) + v_non_cons))));
dt = min(CFL.*(dx./(sqrt(T_non_cons) + v_non_cons)));
%Update old values
rho_old = rho_non_cons;
v_old = v_non_cons;
T_old = T_non_cons;
%Predictor Step
for j = 2:n-1;
dvdx = (v_non_cons(j+1) - v_non_cons(j))/dx;
drhodx = (rho_non_cons(j+1) - rho_non_cons(j))/dx;
dTdx = (T_non_cons(j+1) - T_non_cons(j))/dx;
dlogardx = (log(ar(j+1)) - log(ar(j)))/dx;
%dt is time derivative
%dT is for temprature variable
%Continuity equation
drhodt_p(j) = -rho_non_cons(j)*dvdx - rho_non_cons(j) * v_non_cons(j)*dlogardx - v_non_cons(j)*drhodx ;
%momentum equation
dvdt_p(j) = -v_non_cons(j)*dvdx - (1/gamma)*(dTdx +(T_non_cons(j)/rho_non_cons(j))*drhodx);
%Energy Equation
dTdt_p(j) = -v_non_cons(j)*dTdx - (gamma-1)*T_non_cons(j)*(dvdx + v_non_cons(j)*dlogardx);
%Solution update
rho_non_cons(j) = rho_non_cons(j) + drhodt_p(j)*dt;
v_non_cons(j) = v_non_cons(j) + dvdt_p(j)*dt;
T_non_cons(j) = T_non_cons(j) + dTdt_p(j)*dt;
end
%Corrector Method
for h = 2:n-1;
dvdx = (v_non_cons(h) - v_non_cons(h-1))/dx;
drhodx = (rho_non_cons(h) - rho_non_cons(h-1))/dx;
dTdx = (T_non_cons(h) - T_non_cons(h-1))/dx;
dlogardx = (log(ar(h)) - log(ar(h-1)))/dx;
%Continuity Equation
drhodt_c(h) = -rho_non_cons(h)*dvdx -rho_non_cons(h)*v_non_cons(h)*dlogardx - v_non_cons(h)*drhodx;
%Momentum Equation
dvdt_c(h) = -v_non_cons(h)*dvdx - (1/gamma) * (dTdx+(T_non_cons(h)/rho_non_cons(h))*drhodx);
%Energy Equation
dTdt_c(h) = -v_non_cons(h)*dTdx - (gamma-1)*T_non_cons(h)*(dvdx+v_non_cons(h)*dlogardx);
end
%Computing the Average time derivative
drhodt = 0.5*(drhodt_p+drhodt_c);
dvdt = 0.5*(dvdt_p+dvdt_c);
dTdt = 0.5*(dTdt_p + dTdt_c);
%Final Solution Update
for i = 2:n-1
rho_non_cons(i) = rho_old(i) + drhodt(i)*dt;
v_non_cons(i) = v_old(i) + dvdt(i)*dt;
T_non_cons(i) = T_old(i) + dTdt(i)*dt;
end
%Inlet Boundary Conditions
v_non_cons(1) = 2 * v_non_cons(2) - v_non_cons(3);
%Outlet Boundary Conditions
v_non_cons(n) = 2*v_non_cons(n-1) - v_non_cons(n-2);
rho_non_cons(n) = 2*rho_non_cons(n-1) - rho_non_cons(n-2);
T_non_cons(n) = 2*T_non_cons(n-1) - T_non_cons(n-2);
%Nondimensional mass flow rate
mass_non_cons = rho_non_cons.*ar.*v_non_cons;
%Mach Number
mach_non_cons = v_non_cons./((T_non_cons).^(0.5));
%Nondimensional Pressure
p_non_cons = rho_non_cons.*T_non_cons;
%Throat conditions
%rho_throat2,p_throat2,T_throat2,mach_throat2,v_throat2,mass_throat2
rho_throat2(k) = rho_non_cons(throat_node);
p_throat2(k) = p_non_cons(throat_node);
T_throat2(k) = T_non_cons(throat_node);
mach_throat2(k) = mach_non_cons(throat_node);
v_throat2(k) = v_non_cons(throat_node);
mass_throat2(k) = mass_non_cons(throat_node);
%To compare mass flow rates at different time steps
figure(2)
hold on
if k == 50
plot(x,mass_non_cons,\'b\',\'linewidth\',1.2);
hold on
elseif k == 100
plot(x,mass_non_cons,\'r\',\'linewidth\',1.2);
hold on
elseif k == 150
plot(x,mass_non_cons,\'g\',\'linewidth\',1.2);
hold on
elseif k == 200
plot(x,mass_non_cons,\'y\',\'linewidth\',1.2);
hold on
elseif k == 700
plot(x,mass_non_cons,\'c\',\'linewidth\',1.2);
hold on
elseif k == 1000
plot(x,mass_non_cons,\'.\',\'color\',\'k\',\'markersize\',10);
grid on
title(\'Mass flow rate Non-conservative\');
xlabel(\'Non dimensional distance through nozzle\');
ylabel(\'Non dimensional mass flow\');
legend(\'50dt\',\'100dt\',\'150dt\',\'200dt\',\'700dt\',\'1000dt\');
end
end
end
Outputs
Comparison of time variation at the throat for pressure, mach number, temperature and density along the x axis (length of nozzle) for n = 31
Comparison of Flow variation of pressure, mach number, temperature and density along the x axis (length of nozzle) for n = 31
From the above graph its evident that the flow variation along the length of nozzle for the parameters like pressure, mach number, temperature and density are all same for both Conservative and Non-Conservative flows. Also the mach number at the throat for both Conservative and Non-Conservative flows is 1. Summarising, we can say that all the above mentioned parameters are same for both the cases. Hence both Conservative and Non-Conservative forms give the same results when the flow is Isentropic Supersonic Flow. The only difference is notted in governing equation and the precision.
Comparison of Mass Flow Rate at different time steps for Conservative flow for n = 31
In this graph its evident that the mass flow rate through the nozzle varies very much for timestep 50, in fact it varies more than that of Non-Conservative flow. The reason for this observation is that as the fluid enters the nozzle, there will be a lot of disturbances and hence creating a lot of variations in the flow at initial time stages. As the fluid moves forward in time step the disturbances seems to have decreased (can be observed for timestep 100, 150 and 200). Towards the end of timesteps (timestep 700 & 1000), we can see that the mass flow has reached a steady state.
Comparison of Mass Flow Rate at different time steps for No