Week 7 - Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
AIM:To write a matlab code on Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method.
OBJECTIVE:
- Steady-state distribution of primitive variables inside the nozzle.
- Time-wise variation of the primitive variables.
- Variation of Mass flow rate distribution inside the nozzle at different time steps during the time-marching process.
- Comparison of Normalized mass flow rate distributions of both forms.
EXPLANATION AND THEORY:
Physical Aspects and the Assumptions of the flow:
We consider the steady, isentropic flow through a convergent-divergent nozzle as sketched below. The flow at the inlet to the nozzle comes
from a reservoir where the pressure and temperature are denoted by Po, and To, respectively. The cross-sectional area of the reservoir is large (
theoretically, A -> infinity ), and hence the velocity is very small ( V -> 0). Thus, Po and T0 are stagnation values. There are three important
sections which need to be focussed.
- SubSonic where Mach number is less than 1 in the convergent section part of the nozzle.
- Sonic at the throat (minimum area), where the local velocity of the location is equal to the local speed of sound. (denoted by * values of
P,T,A,V,M). - SuperSonic where the Mach number is greater than 1 The expansion goes isentropically.
Assumptions:
We assume that at a given section, where the cross-sectional area is A, the flow properties are uniform across that section. Hence, although the
area of the nozzle changes as a function of distance along with the nozzle, x, and therefore, in reality, the flow field is two-dimensional (the flow
varies in.the two-dimensional xy space), we make the assumption that the flow properties-vary only with x.Such flow is defined as quasi-onedimensional flow.
Governing Equations(Anlayitical solution):
The governing continuity, momentum, and energy equations for this quasi 1D, steady, isentropic flow can be expressed as:
The governing equations can be derived in two forms:
1) conservative form
2) Non -conservative form.
When we apply the quasi- one dimensional assumptions to these governing equations then we get the equations suitable for present
problem.The governing equation for quasi- one dimensional are
1)Non - Conservative form:
The conservative form of governing equations for quasi 1-D subsonic - supersonic isentropic flow nozzle are:
Continuity equation:
Momentum Equation:
Energy Equation:
Here reservoir temperature and density are denoted by To and p0 , respectively, we define the nondimensional temperature and density, respectively as,
Moreover, L denote the length of the nozzle, we define a dimensionless length as,
Denoting the speed of sound in the reservoir as 'ao' we define a dimensionless velocity as,
Also, the quantity L/ao has the dimension of time, and we define a dimensionless time as,
Finally, we ratio the local area A to the sonic throat area A* and dimensional pressure divided resevoir pressure 'p0' we define a dimensionless
area and presure respectively as,
Then returning to the governing equations and if we introduce these dimensionless variables we get the non-dimensional governing equations
as follows,
Continuity equation:
Momentum Equation:
Energy Equation:
Discretization of equations:
We have finally set up that particular form of the equations that will be most appropriate as well as convenient for the time-marching solution
of quasi-one-dimensional nozzle flow.This method is a simplified form of solving PDE's, as second-order degree equations are hard to solve ,by
the use of this method we get second order accurate time marching without much effort . So, Macormack method is an easy form of explicit
finite-difference technique which is secondorder-accurate in both space and time .It consists of two steps predictor and corrector.
First, consider the predictor step ,In this step we set up the spatial derivatives as forward differences. Also, to reduce the complexity of the
notation, we will drop the use of the prime to denote a dimensionless variable. In what follows, all variables are the nondimensional variables,
denoted earlier by the prime notation.As it is a explicit method ,we find the time derivative terms by using the previous state values and update
the flow field variables .Then we proceed with corrector step.We then calculate the time derivatives by using these updated values by using
backward difference.Then the average value of these time derivatives is found out.
By using this average value we find the corrected values of flow field variables at time t+dt.These steps should be repeated until steady state is
reached:
Calculation of time step (Δt) :
Calculate Δt at all the grid points, i = 1 to i = N, by taking C =0.5, and then choose the minimum value among them.The resulting Δt obtained is
used in the above equations to calculate flow field variables.
Nozzle Shape and Initial Conditions:
The nozzle shape, A = A(x), is specified and held fixed, independent of time, For the case illustrated in this section, we choose a parabolic area
distribution given by,
Therefore, in our choice of initial conditions, we have use any knowledge we may have about a given problem in order to intelligently pick some
initial conditions so that they are choosen close to steady state values so that there will be no chance of solution to go unstable. For example, in
the present problem, we know that density and T decrease and V increases as the flow expands through the nozzle. Hence, we choose initial
conditions that qualitatively behave in the same fashion. For simplicity, let us assume linear variations of the flow-field variables, as a function of
x. For the present case, we assume the following values at time t = 0.
initial conditions at t=0,
ρ = 1-0.3146*x
T =1-0.2314*x
V = (0.1 + 1.09*x)T^1/2
Boundary conditions :
Another aspect of the numerical solution is that of boundary conditions-an all-important aspect, because without the physically proper
implementation of boundary conditions and their numerically proper representation, we have no hope whatsoever in obtaining a proper
numerical solution to our flow problem.
Subsonic inflow boundary (point 1) :
Here, we must allow one variable to float.We choose the velocity V1, because on a physical basis we know the mass flow through the nozzle
must be allowed to adjust to the proper steady state, and allowing V1 to float makes the most sense as part of this adjustment. The value of V1
changes with time and is calculated from information provided by the flow-field solution over the internal points. We use linear extrapolation
from points 2 and 3 to calculate V1 and we take rho(1) =1 and T(1) =1 as fixed values
V(1) = 2*V(2)-V(3)
Supersonic outflow boundary (point N) :
Here, we must allow all flow-field variables to float We again choose to use linear extrapolation based on the flowfield values at the internal
points. Specifically, we have, for the nondimensional variables
Now we write a program using this data and do simulation to find out the flow field variables and required data.
2) Conservative form:
continuity equation:
Momentum Equation:
Energy Equation:
Now converting these equations to non-dimensional form we get the governing equations as,
continuity eqaution:
Momentum Equation:
Energy Equation:
Let us express the equations for quasi-onedimensional flow in generic form. Let us define the elements of the solutions vector U, the flux
vector F, and the source term J as follows.
Now substituting these terms in the above governing equations we get,
To obtain the primitive variables (p, V, T,p, etc.), we must decode the elements U1, U2, and U3 as follows. From the definitions of U1 , U2 , and
U3 given above, we have
In the above equations, we note that the flux vector elements F1, F2, and F3 and source terms are couched in terms of the primitive variables.
When the computer program is written with F 1 , F2 , and F3 expressed directly in terms of p', V', p', and e', instabilities develop during the course
of the time-marching solution .So,we write such that the governing equations are in Pure form in terms of the elements of the solution vector
i.e, in terms of U1,U2 and U3.
Finally the next step to do is discretization .We take the governing equations that consists of flux terms and source terms and discretize them by
using MacCormack's method.
Predictor step: To start the calculation, we use the initial conditions for U1, U2, and U3 to calculate the initial values of F1, F2, and F3 and and weuse these values to calculate the values of U1,U2 and U3 in the next time step .Calculating the time derivatives by using forward difference to
the RHS terms to discretize and these time derivatives are used to calculate the predicted values of U1 ,U2 and U3.
Corrector Step: Now using the predicted values we calculate the new flux terms and source terms .Then we calculate the time derivative terms
by using backward difference.Then we calculate the average values of time derivatives of predictor and corrector steps.
Then we calculate the corrected values of U1,U2,U3 by using average values of time derivative.This is done at every node except the boundary
nodes since we calculate them by using boundary conditions.Then we also calculate primitive variables.These steps are to be repeated until
steady state is obtained.
Initial conditions : The nozzle shape, A = A(x), is specified and held fixed, independent of time, For the case illustrated in this section, we choose a parabolic area distribution.
Since the dependent variables being solved in governing are U1, U2, and U3, we need initial conditions for these same variables at time t = 0 in
order to start the finite-difference solution. The initial conditions for U1, U2, and U3 also allow initial conditions for F1 F2 and F3 to be obtained
from above equations. Such initial conditions for F1, F2, and F3 are needed to form the x derivatives on the right-hand sides of the governing
equations at the first time step.
The initial conditions for U1 , U2, and U3 were synthesized by assuming the following variations for ρ' and T':
These variations slightly more realistic than those assumed in the Nozzle Shape and initial conditions for non conservative form this is in
anticipation that the stability behaviour of the finite difference formulation using the conservation form of the governing equations might be
slightly more sensitive, and therefore it is useful to start with more improve initial conditions than those given for non conservative form. The
initial condition for the variation of V' is synthesized by taking advantage of the fact that one of the dependent variables in our governing
equations, namely U2 is physically the local mass flow; that is, U2= ρ'A'V'. Therefore,for the initial conditions only, let us assume a constant mass
flow through the nozzle and calculate V' as,
The value 0.59 is chosen for U2 because it is close to the exact analytical value of the steady-state mass flow (which for this case is 0.579).
Therefore, the initial condition for V' as a function of x' is obtained by substituting the p' variation given by above equations.Finally, the initial
conditions for U1 , U2, U3 are obtained by substituting the above variations for p', T', and V' into the definitions above
Boundary conditions :
Subsonic inflow boundary (point 1) :
At the subsonic inflow boundary two properties are held fixed and one is allowed to float, and at the supersonic outflow boundary all properties
are allowed to float. In the present formulation, as before, we hold p' and T' fixed at the inflow boundary, both equal to 1.0, and allow V' to float.
By holding p' fixed, then U1 at grid point i = 1 is fixed, independent of time, via U1 = p 'A'. That is,
Supersonic outflow boundary (point N) :
At the supersonic outflow boundary all properties are allowed to float.The flow properties at the downstream, supersonic outflow boundary are
obtained by linear extrapolation from the two adjacent internal points. If N denotes the grid point at the outflow boundary, then:
Now lets write a matlab code which solves the quasi 1D nozzle flow problem using both the conservative and non- conservative forms with the
above data.
MATLAB CODE:
QUASI 1D NOZZLE FLOW CODE:
clear all
close all
clc
%required inputs
L=3;
n=31;
dx=L/(n-1);
C=0.5;
nt=1400;
gamma =1.4;
% solving the governing equations in the non - conservative form
tic
[x,dt,rho,v,T,p,M,flow_rate_n,nt_grid,nn_iter] = non_conservative_form(L,n,C,dx,nt);
toc
time_nc=toc;
fprintf( 'time_of execution for non-conservative form= %0.4gs',time_nc)
% solving governing equations in the conservative form
tic
[dt_c,rho_c,T_c,p_c,M_c,v_c,flow_rate_c,m_analytical,cn_iter,tolerance] = conservative_form(L,n,C,dx,gamma,nt);
toc
time_con=toc;
fprintf( 'time_of execution for conservative form= %0.4gs',time_con)
figure(7)
plot(x, flow_rate_n,'Linewidth',2)
hold on
plot(x, flow_rate_c,'Linewidth',2)
hold on
plot(x,m_analytical,'Linewidth',2)
legend ('non conservative form','conservative form','analytical solution')
xlabel('non dimensional x/L', 'Fontsize', 11);
ylabel('non dimensional Mass flow rate in non conservative form', 'Fontsize',11);
title(sprintf([ 'comparison of the steady-state mass flow variations (on a magnified scale) obtained' ...
'nwith nonconservative and conservation forms of the governing flow equations']))
figure(8)
%bar graph comparing the convergence rate
no_of_iterations = [nn_iter cn_iter];
bar( no_of_iterations ,'g' )
set(gca,'xticklabel', {'non conservative form','conservative form'});
ylabel('no. of iterations');
title(sprintf(['No. of iterations to attain steady state for a tolerance of %g'],tolerance ));
legend('no. of iterations');FUNCTION FOR CONSERVATIVE FORM:
function [dt,rho,T,p,M,v,flow_rate,m_analytical,n_iter,tolerance] = conservative_form(L,n,C,dx,gamma,nt)
% initial conditions
x=linspace(0,L,n);
A=1+2.2*(x-1.5).^2;
throat = (n+1)/2;
nt_grid =linspace(1,nt,nt);
% analytical mass flow rate
m_analytical = 0.579*ones(n);
n_iter =1;
% initial mass flow rate at t=0
U2_i=0.59;
flow_rate_old = 0;
tolerance =1e-4;
error = 999;
for i=1:n
if x(i)>=0 && x(i)<0.5
rho(i)=1.0;
T(i)= 1.0;
else if 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));
else if 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
end
end
%Solution Vectors U1 U2 and U3 which are not primitive variables
for j= 1:n
v(j)=U2_i/(rho(j)*A(j));
U1(j) = rho(j)*A(j);
U2(j) = rho(j)*A(j)*v(j);
U3(j) = rho(j)*A(j)*((T(j)/(gamma-1))+(0.5*gamma*v(j)^2));
end
% calculating minimum of dt calculated at all nodes
for m=1:n
dt_x(m)= C*(dx/(T(m)^(1/2)+v(m)));
end
% time step
dt=min(dt_x);
U2_initial= U2;
for k = 1:nt
U1_old = U1;
U2_old = U2;
U3_old = U3;
%flux terms which are functions of solution vectors
F1 = U2;
F2 = ((U2.^2)./(U1))+ (((gamma-1)/gamma).*(U3 - ((gamma/2).*(((U2.^2)./U1)))));
F3 = ((gamma.*U2.*U3)./U1) - (gamma.*(gamma-1)/2)*((U2.^3)./U1.^2);
for i = 2:n-1
dF1_dx = (F1(i+1) - F1(i))/dx;
dF2_dx = (F2(i+1) - F2(i))/dx;
dF3_dx = (F3(i+1) - F3(i))/dx;
%Continuity equation predictor step
dU1_dt_p(i) = -dF1_dx;
%Momentum Equation
J2(i)= ((1/gamma)*rho(i)* T(i)*((A(i+1)-A(i))/dx));
dU2_dt_p(i) = -dF2_dx + J2(i);
%energy equation of predictor step
dU3_dt_p(i) = -dF3_dx;
%%Solution Update of solution vectors
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
%updation of initial values
rho = U1./A;
v = U2./U1;
T = (gamma-1).*((U3./U1)-(0.5*gamma.*v.^2));
p = rho.*T;
% computing new set of values using predicted values
F1 = U2;
F2 = ((U2.^2)./(U1))+ (((gamma-1)/gamma).*(U3 - ((gamma/2).*(((U2.^2)./U1)))));
F3 = ((gamma.*U2.*U3)./U1) - (gamma.*(gamma-1)/2).*((U2.^3)./U1.^2);
% corrector step
% calculating derivative terms using backward difference for corrector method
for j = 2:n-1
dF1_dx = (F1(j) - F1(j-1))/dx;
dF2_dx = (F2(j) - F2(j-1))/dx;
dF3_dx = (F3(j) - F3(j-1))/dx;
%Continuity equation corrector step
dU1_dt_c(j) = -dF1_dx;
%Momentum Equation
J2(j)= ((1/gamma)*rho(j)* T(j)*((A(j)-A(j-1))/dx));
dU2_dt_c(j) = -dF2_dx + J2(j);
%energy equation of corrector step
dU3_dt_c(j) = -dF3_dx;
end
%average values
dU1_dt_avg = 0.5*(dU1_dt_p+dU1_dt_c);
dU2_dt_avg = 0.5*(dU2_dt_p+dU2_dt_c);
dU3_dt_avg = 0.5*(dU3_dt_p+dU3_dt_c);
%final corrected Values
for m = 2:n-1
U1(m) = U1_old(m)+ dU1_dt_avg(m)*dt;
U2(m) = U2_old(m)+ dU2_dt_avg(m)*dt;
U3(m) = U3_old(m)+ dU3_dt_avg(m)*dt;
end
% inlet boundary conditions
U1(1)= rho(1)*A(1);
U2(1) = 2*U2(2)-U2(3);
v(1) = U2(1)/(rho(1)*A(1));
U3(1)= U1(1)*((T(1)/(gamma-1))+((gamma*0.5)*v(1)^2));
% Outlet conditions with all the three float
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);
% update of primitive variables
rho = U1./A;
v= U2./U1;
T = (gamma-1).*((U3./U1)-(0.5*gamma.*v.^2));
% calculating pressure(p), mach number (M)and mass flow rate
p = rho.*T;
M= (v./((T).^0.5));
flow_rate = rho.*v.*A;
%capturing flow field values at throat for each time step
rho_throat(k) = rho(throat);
v_throat(k) = v(throat);
T_throat(k) = T(throat);
p_throat(k) = p(throat);
M_throat(k) = M(throat);
flow_rate_throat(k) = flow_rate(throat);
error = max(abs(flow_rate-flow_rate_old));
if error > tolerance
n_iter = n_iter+1;
flow_rate_old =flow_rate;
end
% non dimensional mass flow distribution across the nozzle at different time steps
% compared with analytical mass flow rate solution
figure(4)
if k == 50
plot(x, flow_rate,'Linewidth', 2);
hold on
else if k == 100
plot(x, flow_rate, 'Linewidth',2);
hold on
else if k == 200
plot(x, flow_rate, 'Linewidth',2);
hold on
else if k == 500
plot(x, flow_rate, 'Linewidth',2);
hold on
else if k == 700
plot(x, flow_rate,'Linewidth',2);
hold on
plot(x, m_analytical,':')
hold on
plot(x, U2_initial,'--')
legend('50Deltat','100Deltat','200Deltat','400Deltat','700Deltat','analytical solution','initial flow rate');
title('non dimensional mass flow distribution across the nozzle at different time steps before steady state-conservative form');
xlabel('non dimentional x/L', 'Fontsize', 11);
ylabel('Mass flow rate conservative', 'Fontsize',11);
axis([0 3 0.5 1])
end
end
end
end
end
end
% plotting steady state values of flow field variables at each node
figure(5)
sgtitle(sprintf('steady state values of flow field variables in conservative form ,number of timestep=%d',k))
subplot (5,1,1)
plot(x,M,'b','Linewidth',2);
xlabel('non dimensional length');
ylabel('M');
legend('mach number')
subplot (5,1,2)
plot(x,v,'r','Linewidth',2);
xlabel('non dimensional length');
ylabel('V/a_o');
legend('non dimensional velocity')
subplot (5,1,3)
plot(x,rho,'b','Linewidth',2);
xlabel('non dimensional length');
ylabel('ρ/ρ_o','FontSize', 10);
legend('non dimensional density')
subplot (5,1,4)
plot(x,T,'y','Linewidth',2);
xlabel('non dimensional length','FontSize', 10);
ylabel('T/T_o','FontSize', 10);
legend('non dimensional temperature')
subplot (5,1,5)
plot(x,p,'b','Linewidth',2);
xlabel('non dimensional length');
ylabel('p/p_o');
legend('non dimensional pressure')
% timewise variation of the flow field variables and mach number at the nozzle throat
figure(6)
plot (nt_grid,M_throat,'Linewidth',2)
hold on
plot (nt_grid,rho_throat, 'Linewidth' ,2)
hold on
plot(nt_grid,v_throat,'Linewidth',2)
hold on
plot (nt_grid,T_throat,'Linewidth',2)
hold on
plot (nt_grid,p_throat,'Linewidth',2)
legend('mach number','non dimensional density(ρ/ρ_o)', 'non dimensional velocity(V/a_o)'...
,'non dimensional temperature(T/T_o)','non dimensional pressure(p/p_o)')
xlabel('Number of time steps')
axis([0 1400 0.5 1.2])
title('timewise variation of non dimensional flow field variables at throat in conservative form')
endFUNCTION FOR NON - CONSERVATIVE FORM:
function [x,dt,rho,v,T,p,M,flow_rate_n,nt_grid,n_iter] = non_conservative_form(L,n,C,dx,nt)
% initial conditions
x=linspace(0,L,n);
rho =1-0.3146*x;
T= 1-0.2314*x;
v=(0.1+1.09*x).*T.^0.5;
A=1+2.2*(x-1.5).^2;
gamma =1.4;
throat = (n+1)/2;
nt_grid =linspace(1,nt,nt);
initial_flow =rho.*v.*A;
flow_rate_old = 0;
tolerance =1e-4;
error = 999;
n_iter=1;
% analytical mass flow rate
m_analytical = 0.579*ones(n);
% calculating minimum of dt calculated at all nodes
for m=1:n
dt_x(m)= C*(dx/(T(m)^(1/2)+v(m)));
end
% time step
dt=min(dt_x);
for k= 1:nt
rho_old = rho;
T_old = T;
v_old = v;
% predictor method
% calculating derivative terms using forward difference for predictor method
for i= 2:n-1
dvdx =(v(i+1)-v(i))/dx;
drhodx = (rho(i+1)-rho(i))/dx;
dT_dx=( T(i+1)-T(i))/dx;
dlogA_dx = (log(A(i+1))-log(A(i)))/dx;
%continuity equation
drhodt_p(i) = -rho(i)*dvdx-rho(i)*v(i)*dlogA_dx-v(i)*drhodx;
%momentum equation
dvdt_p(i) =-v(i)*dvdx- (1/gamma)*(dT_dx+(T(i)/rho(i))*drhodx);
%energy equation
dTdt_p(i) = -v(i)*dT_dx - (gamma-1)*T(i)*(dvdx + (v(i)*dlogA_dx));
% predictor solution step update
rho(i)= rho(i)+drhodt_p(i)*dt;
v(i)= v(i)+dvdt_p(i)*dt;
T(i)= T(i)+dTdt_p(i)*dt;
end
% corrector step
% calculating derivative terms using backward difference for corrector method
for i= 2:n-1
dvdx =(v(i)-v(i-1))/dx;
drhodx = (rho(i)-rho(i-1))/dx;
dT_dx=( T(i)-T(i-1))/dx;
dlogA_dx = (log(A(i))-log(A(i-1)))/dx;
%continuity equation
drhodt_c(i) = -rho(i)*dvdx-rho(i)*v(i)*dlogA_dx-v(i)*drhodx;
%momentum equation
dvdt_c(i) =-v(i)*dvdx- (1/gamma)*(dT_dx+(T(i)/rho(i))*drhodx);
%energy equation
dTdt_c(i) = -v(i)*dT_dx - (gamma-1)*T(i)*(dvdx + (v(i)*dlogA_dx));
end
% compute the average time derivatives
dvdt = 0.5*(dvdt_p+dvdt_c);
dTdt=0.5*(dTdt_p+dTdt_c);
drhodt =0.5*(drhodt_p+drhodt_c);
% final step update
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
% Inlet Boundary Conditions
v(1) = 2*v(2) - v(3);
% Outlet Boundary Conditions
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);
% calculating pressure(p), mach number (M)and mass flow rate
p=rho.*T;
M = (v./((T).^0.5));
flow_rate_n = rho.*v.*A;
%capturing flow field values at throat for each time step
rho_throat(k) = rho(throat);
v_throat(k) = v(throat);
T_throat(k) = T(throat);
p_throat(k) = p(throat);
M_throat(k) = M(throat);
flow_rate_throat(k) = flow_rate_n(throat);
error = max(abs(flow_rate_n-flow_rate_old));
if error > tolerance
n_iter = n_iter+1;
flow_rate_old =flow_rate_n;
end
figure(1)
% non dimensional mass flow distribution across the nozzle at different time steps
% compared with analytical mass flow rate solution
if k == 50
plot(x, flow_rate_n,'Linewidth', 2);
hold on
else if k == 100
plot(x, flow_rate_n, 'Linewidth',2);
hold on
else if k == 200
plot(x, flow_rate_n, 'Linewidth',2);
hold on
else if k == 500
plot(x, flow_rate_n, 'Linewidth',2);
hold on
else if k == 700
plot(x, flow_rate_n,'Linewidth',2);
hold on
plot(x, m_analytical,':')
hold on
plot(x, initial_flow,'--')
legend('50Deltat','100Deltat','200Deltat','400Deltat','700Deltat','analytical solution','initial flow rate');
title('non dimensional mass flow distribution across the nozzle at different time steps before steady state- non conservative form');
xlabel('non dimensional x/L', 'Fontsize', 11);
ylabel('Mass flow rate in non conservative form', 'Fontsize',11);
end
end
end
end
end
end
% plotting steady state values of flow field variables at each node
figure(2)
sgtitle(sprintf('steady state values of flow field variables in non conservative form , number of time steps=%d',k))
subplot (5,1,1)
plot(x,M,'b','Linewidth',2);
xlabel('non dimensional length');
ylabel('M');
legend('mach number')
subplot (5,1,2)
plot(x,v,'r','Linewidth',2);
xlabel('non dimensional length');
ylabel('V/a_o');
legend('non dimensional velocity')
subplot (5,1,3)
plot(x,rho,'b','Linewidth',2);
xlabel('non dimensional length');
ylabel('ρ/ρ_o','FontSize', 10);
legend('non dimensional density')
subplot (5,1,4)
plot(x,T,'y','Linewidth',2);
xlabel('non dimensional length','FontSize', 10);
ylabel('T/T_o','FontSize', 10);
legend('non dimensional temperature')
subplot (5,1,5)
plot(x,p,'b','Linewidth',2);
xlabel('non dimensional length');
ylabel('p/p_o');
legend('non dimensional pressure')
figure(3)
% timewise variation of the flow field variables and mach number at the nozzle throat
plot (nt_grid,M_throat,'Linewidth',2)
hold on
plot (nt_grid,rho_throat,'Linewidth',2)
hold on
plot(nt_grid,v_throat,'Linewidth',2)
hold on
plot (nt_grid,T_throat,'Linewidth',2)
hold on
plot (nt_grid,p_throat,'Linewidth',2)
legend('mach number','non dimensional density(ρ/ρ_o)', 'non dimensional velocity(V/a_o)'...
,'non dimensional temperature(T/T_o)','non dimensional pressure(p/p_o)')
xlabel('Number of time steps')
title('timewise variation of non dimensional flow field variables at throat in non conservative form')
endRESULTS:
1) Steady state values of flow field variables in non conservative form,
2) Steady state values of flow field variables in conservative form,
3) In the above plots you case see the timewise variation of non dimensional flow field variables at throat in both non conservative form and conservative form respectively attaining steady state.You can also observe that the oscillations reduced almost ,means attaining steady state ,in non conservative form you can see that ends almost nearly after 450 timesteps but for the conservative form the oscillations continued upto 600 time steps but after that it attained steady state. By this you can note that non conservative form reaches steady state earlier than conservative form.
4) In the above plots you can observe the non dimensional mass flow distribution across the nozzle at different time steps in both non conservative form and conservative form respectively,
By observing the above plot you can conclude when the mass is conserved,in both the forms you can observe that it is conserved almost after 700 times steps
5) In the below result we have compared the normalized mass flow rate i.e, steady state mass flow variations obtained with non conservative and conservative forms of governing equations with the analytical solution across the nozzle.
The comparison shown in above figure illustrates a general advantage of the conservation form of the equations. 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,the mass flow is a primary result from these equations. In contrast, the dependent variables in the non conservation form of the equations are the primitive variables, and the mass flow is obtained only as a secondary result.
6) In the below plot we have compared between both the forms showing that which form is attaining steady state faster to attain a tolerance of 1e-4.
Therefore by observing the above bar graph you can conclude that non conservative form attains steady state faster than conservative form and you can stop the simulation when the steady state is reached.So finally you can conclude that non conservative form is faster than conservative and you can time of execution in below figure:
GRID DEPENDANCY TEST:
Let us take a look at the primitive variables,rho, temperature, pressure, and Mach number at the nozzle throat, as tabulated in below table:
The first row gives the exact, analytical results. The second and third rows give the numerical results for the nonconservation and the conservation forms at n=31, respectively.Note that the nonconservation form results are distinctly closer to the exact values. The The fourth and fifth rows give the numerical results for the nonconservation and the conservation forms at n=61, respectively,for a grid with twice as many grid points (61 in comparison to 31 points). A comparison of the rows with n=31 to n=61 of both the forms are an indication of grid independence for the conservation form and non-conservation form results since the change in the values are very small.
CONCLUSION:
1) The conservation form yields a better mass flow distribution, since the mass flow itself is a dependent value.
2) The nonconservation form seems to produce slightly more accurate results for the primitive variables, and the conservation form seems to produce slightly more accurate results for the flux variables.
3) Non-conservative form attains steady state faster than the conservative form.
4) Both the forms conservative and non-conservative are grid independent.