Week 7 - Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
Objective:
1)Solving Quasi 1D Nozzle Flows through Conservative and Non-Conservative Methods.
2)Finding values of flow field variables across nozzle after attaining steady state for both the forms.
3)To capture the timewise varation of flow field variables at the nozzle throat.
4)Distributions of the nondimensional mass flow rate as a function of distance through the nozzle at different time steps compared with the analytical solution for both the forms.
5)Comparison of the steady-state mass flow variations obtained with the nonconservation and conservation forms of the governing flow equations with the analytical solution.
6)To perform grid dependence test and find the fastest of the two methods which attain steady state first and also to find which method executes the results faster.
(Theory - John D Anderson COMPUTATIONAL FLUID DYNAMICS)
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
1) SubSonic where Mach number is less than 1 in the convergent section part of the nozzle.
2) 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)
3) 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-one-dimensional flow.
Governing Equations(Anlayitical solution): The governing continuity, momentum, and energy equations for this quasi 1D, steady, isentropic flow can be expressed, respectively, as
Solving them analytical we arrive at the below graphs.
What we want to accomplish here is to illustrate the application of various CFD techniques, and we are intentionally choosing a flow problem with a known analytic solution for this illustration.
CFD Solution :
the setup :
In this section, we will set up three eschelons of equations as follows
1. The governing flow equations will be couched in terms of partial differential equations suitable for the time-marching solution of quasi-one-dimensional flow (the closed-form algebraic equations discussed above for a steady flow and are not suitable for the present purpose).
2. The finite-difference expressions pertaining to MacCormack's technique as applied to this problem will be set up.
3. Other details for the numerical solution such as the calculation of the time step and the treatment of boundary conditions will be formulated.
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:
At this stage, we could readily proceed to set up our numerical solution of above equations . Note that these are written in terms of dimensional variables.This is fine, and many CFD solutions are carried out directly in terms of such dimensional variables. Indeed, this has an added engineering advantage because it gives you a feeling for the magnitudes of the real physical quantities as the solution progresses. However, for nozzle flows, the flow-field variables are frequently expressed in terms of non dimensional variables, where the flow variables are referenced to their reservoir values,(where the dimensional flow variables are divided by their respective reservoir values). The nondimensional variables p/p0 , p/p0 , and T/T0 vary between 0 and 1, which is an "aesthetic" advantage when presenting the results. Because fluid dynamicists dealing with nozzle flows so frequently use these nondimensional terms, we will follow suit here.
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 second-order-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 equation :
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 we use 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 given by
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.
% Main program
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 to solve non conservative form of governing equations
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('50\Deltat','100\Deltat','200\Deltat','400\Deltat','700\Deltat','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')
end
% Function to solve conservative form of governing equations
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('50\Deltat','100\Deltat','200\Deltat','400\Deltat','700\Deltat','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 time step=%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')
endResults :
1) Steady state values of flow field variables in non conservative form,
2) Steady state values of flow field variables in conservative form,
The above plots are almost same , when compared with the analytical solution from John D Anderson book.
3) In the below 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 below 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 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.
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.
Grid Dependency 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. We can also note that, by doubling the number of grid points, the steady-state numerical results are slightly closer to the exact, analytical values .You can also observe that even after having 61 grid points for conservative form still the results are not as close as the nonconservation form results with half as many grid points. We can finally conclude that the solution from both the forms are grid independent.
Conclusion :
After the completion of entire simulation and seeing the results we can finally conclude that in summary, for the given flow problem, we cannot establish a clear superiority of the conservation form of the governing equations over the nonconservation form. In essence, from all our previous discussions, we can only make the following observations:
1. The conservation form yields a better mass flow distribution. The conservation form simply does a better job of conserving mass since the mass flow itself is a dependent value.
2.There is no clear superiority of either form in terms of accuracy of the results. 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. The results in either case are certainly satisfactory.
3 .Comparing the amount of calculational effort to achieve a solution, we note that the solution of the conservation form requires marginally more work. Most of this is due to the need to decode the primitive variables from the flux variables. such decoding is not necessary when you are solving the nonconservation form.
4.Non-conservative form attains steady state faster than the conservative form.
5.Both the forms conservative and non-conservative are grid independent.