Week 7 - Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
Problem specification:
In this challenge, you are going to simulate the isentropic flow through a quasi 1D subsonic-supersonic nozzle. You will derive both the conservation and non-conservation forms of the governing equations and solve them using the MacCormack's technique. You need to determine the steady-state temperature distribution for the flow-field variables and investigate the difference between the two forms of governing equations by comparing their solutions.
Need to show the following plots inside the report
1. Steady-state distribution of primitive variables inside the nozzle
2. Time-wise variation of the primitive variables
3. Variation of Mass flow rate distribution inside the nozzle at different time steps during the time-marching process
4. Comparison of Normalized mass flow rate distributions of both forms
Expected Results:
- Perform a grid independence test on the solution and explain all your findings.
- The code must be modular. You need to write separate functions to solve both forms.
- Discuss all the plots you are providing with observed inferences.
THEORY:
MacCormack technique:
It is an explicit finite difference technique which is second-order accurate in both space and time.
Let's consider an unsteady,2-d isentropic flow
Eulers equations in non-conservative form:
Continuity: `(deltarho)/(
x-momentum:`(deltau)/(deltat)=
y-momentum:`(deltav)/(deltat)=
Energy:`(deltae)/(deltat)=-(u(
In the above equation, we assumed nobody forces and no volumetric heat addition i.e f=0 and q=0 the above equation is hyperbolic w.r.t time we assume that flow field at the time (t) is known and we will identify at (t+Deltat).
Here we are choosing density rho as our dependent flow variable for which we are going to derive the MacCormack method,
Let `(rho^(t))_(i,j)` denote the density at grid points (i,j) at the time(t) density at the same grid point at the time (t+Deltat) is given by,
`(rho^(t+Deltat))_(i,j)=(rho^(
Where `((deltarho)/(deltat))_(av` denotes the mean value of `(deltarho)/(deltat)` between the times(t) and (t+Deltat)
values of `((deltarho)/(deltat))_(av` is calculated to preserve second-order accuracy without the need to calculate values of the second-order time derivative which is the term involves a lot of algebra.
Similar relations are written for the other flow fiel variables.
`(u^(t+Deltat))_(i,j)=(u^(t))_
`(v^(t+Deltat))_(i,j)=(v^(t))_
`(e^(t+Deltat))_(i,j)=(e^(t))_
From equation (5);
Average time derivative `(deltarho)/(
Predictor step:
In the continuity equation (1) discretize the spatial derivatives on the right-hand side using forward difference scheme.
Predicted value of density (rho)t+Deltat from the first two terns of taylor series as follows,
`(rho^(t+Deltat))_(i,j)=(rho^(
similarly:
`(u^(t+Deltat))_(i,j)=(u^(t))_
`(v^(t+Deltat))_(i,j)=(v^(t))_
`(e^(t+Deltat))_(i,j)=(e^(t))_
In the above equations, 10 to 13 numbers for the time derivative on the right-hand side are obtained from the equation.
Corrector step:
In the corrector step,we use the calculated predicted value of the time derivative at the time `(t+Deltat).((deltarho)/(
`(deltarho)/(deltat)=-(rho^-t+
The average value of the time derivative of density in equation (5) is obtianed from the arithmetic mean of `((deltarho)/(deltat))^t)_(
this allows us to obtain the final corrected valued of density at the time (t+deltat) from equation 5
`(rho^(t+Deltat))_(i,j)=(rho^(
A similar procedure can be applied to calculate momentum and energy.
Let `(rho^(t))_(i,j)` denote the density at grid points (i,j) at the time(t) density at the same grid point at the time (t+Deltat) is given by,
`(rho^(t+Deltat))_(i,j)=(rho^(
Where `((deltarho)/(deltat))_(av)` denotes the mean value of `(deltarho)/(deltat)` between the times(t) and (t+Deltat)
values of `((deltarho)/(deltat))_(av)` is calculated to preserve second-order accuracy without the need to calculate values of the second-order time derivative which is the term involves a lot of algebra.
Similar relations are written for the other flow fiel variables.
`(u^(t+Deltat))_(i,j)=(u^(t))_
`(v^(t+Deltat))_(i,j)=(v^(t))_
`(e^(t+Deltat))_(i,j)=(e^(t))_
From equation (5);
Average time derivative `(deltarho)/(
Predictor step:
In the continuity equation (1) discretize the spatial derivatives on the right-hand side using forward difference scheme.
Predicted value of density (rho)t+DEltat from the first two terns of taylor series as follows,
`(rho^(t+Deltat))_(i,j)=(rho^(
similarly:
`(u^(t+Deltat))_(i,j)=(u^(t))_
`(v^(t+Deltat))_(i,j)=(v^(t))_
`(e^(t+Deltat))_(i,j)=(e^(t))_
In the above equations, 10 to 13 numbers for the time derivative on the right-hand side are obtained from the equation.
Corrector step:
In the corrector step,we use the calculated predicted value of the time derivative at the time `(t+Deltat).((deltarho)/(
`(deltarho)/(deltat)=-(rho^-t+
The average value of the time derivative of density in equation (5) is obtianed from the arithmetic mean of `((deltarho)/(deltat))^t)_(
this allows us to obtain the final corrected valued of density at the time (t+deltat) from equation 5
`(rho^(t+Deltat))_(i,j)=(rho^(
A similar procedure can be applied to calculate momentum and energy.
Numerical solution of quasi one-dimensional supersonic nozzle flow using MacCormack method:
Assumptions: the flow at the inlet to the nozzle comes from a reservoir where the pressure, temperature and density are denoted by the cross-sectional area of the reservoir is large (A.infinity)and hence velocity is very small (V>0).
We consider the steady, isentropic flow through a convergent-divergent nozzle as sketched in The flow at the inlet to the nozzle comes from a reservoir where the pressure and temperature are denoted by p0 , and To, respectively. Across-sectional area of the reservoir is large ( theoretically, A >infinity ), and hence the velocity is very small ( V> O). Thus, Po and T0 are the stagnation values,_ or total
pressure and total temperature, respectively. The flow expands 1sentrop1cally. to supersonic speeds at the nozzle exit, where the exit pressure, temperature, velocity and Mach number are denoted by Pe, Te, V.,, and Me, respectively. The flow IS locally subsonic in the convergent section if the nozzle, some at the throat (minimum area), and supersonic at the d1vergent section. flow
(M = 1) at the throat means that the local velocity at this location 1s equal to the local speed of sound.
Using an asterisk to denote sonic flow values, we have at the throat V = V* = a*. Similarly, the sonic flow values of pressure and temperature are denoted by p* and T*, respectively. The area of the sonic throat is denoted by A*. 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 .ry space), we make the assumption that the flow properties-vary only with x; this is tantamount to assuming uniform flow properties across any given cross-section. Such flow is defined as quasi-one-dimensional flow. The governing continuity, momentum, and energy equations for this quasi one-dimensional, steady, isentropic flow can be expressed, respectively, as
Continuity:
Conservative form: `(delta)/(deltat)(rhoA)+
Non conservative form:(delta)/(deltat)(rhoA)+V(
Momentum:
Conservtaive:`(delta)/(deltat)
Non conservative form:`rho(deltaV)/(deltat)+
Energy:
conservative:`(delta)/(deltat)(
Non conservative:`rho(deltae)/(
The reason for obtaining the energy equation in the form of equation 21 is that for a calorically perfect gas it leads directly to form of the energy equation in terms of temperature T.for our solution of the quasi one dimensional nozzle flow of a calorically perfect gas this is a fundamental variable T and therefore it is convenient to deal with it as the primary dependent variable in the energy section.
For a calorically perfect gas e= CvT
hence the energy equation of conservative and non conservative form from equation 20 and 21 becomes
`conservative form:C_V(delta)/(deltat)(
`Non conservative:rhoC_v(deltaT)/(
Momentum:
`conservtaive:(delta)/(deltat)
`non conservative form:rho(deltaV)/(deltat)+
Non-dimensional variables:
Dimensionless temperature T' =T/To
Dimensionless density rho'=rho/rho0
where To and rho0 represent reservoir temperature and density respectively.
dimensionless length x'=x/L
where L represents the length of the nozzle.
speed of sound in the reservoir ao=squareroot of grammar to.
dimensionless velocity V'=V/a0
dimensionless time t'=T/L/ao.
dimensionless area A'=A/A local area /sonic throat area.
Returning to our derived continuity momentum and energy equations and substituting non-dimensional variables.
Continuity:
conservative form:`(delta)/(deltat')(rho'A'
non-conservative:`(deltarho')/
Momentum:
Conservative form:`(delta)/(deltat')(rho'A'
Non conservative form:`(deltaV')/(deltat)=-V'(
Energy:
conservative:`(delta)/(deltat'
Non conservative:`(deltaT')/(
Meshing:
The first grid point labelled point 1 is assumed to be in the reservoir points are evenly distributed along the x-axis with Deltax grid spacing point(i) is simply an arbitrary grid point with points(i-1) and (i+1) as the adjacent points
For this particular problem, we are taking the number of grid points to be 31 and 61(as to run the grid dependence test)
Average time derivatives:
`(rho^(t+Deltat))_(i,)=(rho^(
`(V^(t+Deltat))_(i,)=(V^(t))_(
`(T^(t+Deltat))_(i,)=(T^(t))_(
Initial conditions(at the time =0);
Depending on the values of x-node we are setting up the values for rho,V and T.
rho=1-0.3146x.
V=1-0.2314x
T=(0.1=1.09x)T^0.5.
Nozzle shape:
For area(A), we are going to set up a profile.
A=1+2.2(x-1.5)^2 0<x<3.
Mach number:
`M=v/sqrt(t)`
Calculation of time step:
For both non-conservative and conservative form of governing equation, we are going to take the time step.
`Deltat=c(Deltax)/(a+v)`-----
where C=courant number
V=local flow velocity
a=locl speed of sound.
for an explicit numerical solution to be stable C<1, here we are taking C=0.5.
Equation 40 is the CFL criterion for a one-dimensional flow along with C<1 it simply states that delta t must be less than or at best equal to the time taken by a sound wave to move from one grid point to the next.
Deltax is the same throughout the flow, but both V' and a are variables. hence at a given grid points at a given time step.
`(Deltat^t)_(i)=C(Deltax)/((a^
And at an adjacent grid point:
`(Deltat^t)_(i+1)=C(Deltax)/((
Boundary conditions:
subsonic inflow boundary (at node1 )
We must allow one variable to float, we choose the velocity V because on a physical basis we know that the mass flow through the nozzle must be allowed to adjust to the proper steady-state and allowing V to float makes the complete sense as part of this adjustment
V1=2V2-V3
Density and temperature are fixed (independent of time)
rho1=1 and T1=1
Supersonic outflow boundary(at node N)
we must allow all flow field variables to float.
Vn=2Vn-1-Vn-2
rhon=2rhon-1-rhon-2
Tn=2Tn-1-Tn-2
Conservative form of governing equations:
Governing equations for quasi-one-dimensional flow in conservative form is given by
Continuity:
conservtaive form;`(delta)/(deltat')(rho'A'
Momentum:
conservative form:`(delta)/(deltat')(rho'A'
Energy:
conservtaive form:`(delta)/(deltat')[rho'(
The governing equations for quasi one-dimensional flow can be expressed in a more generic form by defining the elements as the solution vector U the flux vector F and the source term J as follows.
U1=rho'A'
U2=rho'A'V'
U3=`rho'[(e')/(gamma-1)+(
F1=rho'A'V'
F2=rho'A'(V')^2+1/gamma*rho'A'
F3=`rho'(e')/(gamma-1)+(gamma)
J2=`1/gammaP'(deltaA')/(
with these elements governing equations 51 and 53 it can be written as;
Continuity:
`(deltaU1)/(deltat')=-(
Momentum:
`(deltaU2)/(deltat')=-(
Energy:
`(deltaU3)/(deltat')=-(
In the conservative form of governing equations the dependent variables to obtain the primitive variables(rho,V,P,T etc...) we must decode the elements U1,U2,U3 as follows.
rho'=U1/A'------57
V'=U2/U1-----58
T'=`e'-(gamma-1)(U3)/(U2)-(
P'=rho'T'-------60
Therefore after we obtain U1, U2 and U3 at each time step from the numerical solution of equations 54 to 56 we can immediately calculate the corresponding primitive variables at each time step rho', V', T' and P' from equations 57 to 60.
Pure form of flux vector is given by;
F1=U2
F2=`(U2^2)/(U1)+(gamma-1)/(
F3=`gamma(U2U3)/(U1)-gamma(y-
J2=`(gamma-1)/(gamma)(U3-(
This is the pure form of governing equation in conservative form.
Initial conditions:
for 0<x'0.5: rho'=1.0 and T' =1.0
for 0.5<1.5: rho'=1.0-0.366(x'-0.5) and T'=1.0-0.167(x'-0.5)
for 1.5<x'<3.0: rho'=0.634-0.3879(x'-1.5) and T'=0.833-0.3507(x'-1.5) .
for the initial condition only,lets assume a constant mass flow through the nozzle and calculate V as,
`V'=(U2)/(rho'A')=(0.59)/(rho'
The value of 0.59 is chosen for U2 because it is the closest to the exact analytical value of the steady state mass flow (i.e)=0.579.
Boundary conditions:
Subsonic inflow boundary (at node 1)
Two properties are held fixed and one is allowed to float rho and T are fixed at inflow boundary both equal to 1 and allowing V to float.
`(U1_1)_(i=1)=rho'A'_(i=1)`
`(U2)_(i=2)=2(U2)_(i=2)-(U2)_(
`(U3)_(i=1)=rho'[(e')/(gamma-
Supersonic outflow boundary(at node N)
All properties are allowed to float:
`(U1)_n=2(U1)_(n-1)-(U1)_(n-1)`
`(U2)_n=2(U2)_(n-1)-(U2)_(n-1)`
`(U3)_n=2(U3)_(n-1)-(U3)_(n-1)`
Main code :
%Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
close all
clear all
clc
%meshing
n=31; %no. of grid points
x=linspace(0,3,n); %grid points distribution
dx=x(2)-x(1); %spacing between the grid points
gamma=1.4;
a=1+2.2*(x-1.5).^2; %area
Throat=find(a==1); %finding the throat
nt=5400; %no. of time steps
%function for the Non conservation form
[mass_flow_rate_1] = non_conservativr_1D(a,n,x,dx,nt,gamma,Throat);
%function for the conservation form
[mass_flow_rate_2]=conservative_1D(a,n,x,dx,nt,gamma,Throat);
figure('Name','comparsion between conservation,non_conservation and exact analytical solution')
plot(x,mass_flow_rate_1,'k')
hold on
plot(x,mass_flow_rate_2,'r')
hold on
plot(x,0.529,'p')
grid on
xlabel('length of the nozzle')
ylabel('mass flow rate')
legend('Non conservation form','Conservation form','exact analytical solution ')
Calling the function for non conservation form:
function [mass_flow_rate_1]=non_conservativr_1D(a,n,x,dx,nt,gamma,Throat)
%to calculate initial profiles
rho=1-0.3146*x; %density
t=1-0.2314*x; %temperature
v=(0.1+1.09*x).*t.^0.5; %velocity
dt=0.0226;
%outer time loop
for k=1:nt
rho_old=rho;
v_old=v;
t_old=t;
%predictor method
for j=2:n-1
dvdx=(v(j+1)-v(j))/(dx);
drhodx=(rho(j+1)-rho(j))/(dx);
dtdx=(t(j+1)-t(j))/dx;
dlogadx=(log(a(j+1))-log(a(j)))/(dx);
%continuity equation
drhodt_p(j)=-rho(j)*dvdx-rho(j)*v(j)*dlogadx-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)*dlogadx);
%solution update
rho(j)=rho(j)+drhodt_p(j)*dt;
v(j)=v(j)+dvdt_p(j)*dt;
t(j)=t(j)+dtdt_p(j)*dt;
end
%corrector method
for s=2:n-1
dvdx=(v(s)-v(s-1))/(dx);
drhodx=(rho(s)-rho(s-1))/(dx);
dtdx=(t(s)-t(s-1))/dx;
dlogadx=(log(a(s))-log(a(s-1)))/(dx);
%continuity equation
drhodt_c(s)=-rho(s)*dvdx-rho(s)*v(s)*dlogadx-v(s)*drhodx;
%momentum equation
dvdt_c(s)=-v(s)*dvdx-(1/gamma)*(dtdx+(t(s)/rho(s))*drhodx);
%energy equation
dtdt_c(s)=-v(s)*dtdx-(gamma-1)*t(s)*(dvdx+v(s)*dlogadx);
end
%computee 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 value 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
%apply the boundary condition
% AT subsonic iNLET flow
v(1)=2*v(2)-v(3);
%at supersonic outflow
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);
%static pressure
P=rho_old.*t_old;
% mach no.
M=v_old./(t_old.^0.5);
% mass flow rate
mass_flow_rate_1=rho_old.*a.*v_old;
% Pressure, temperature, mach number,at the nozzle throat
P_th(k)=P(Throat);
rho_th(k)=rho(Throat);
M_th(k)=M(Throat);
t_th(k)=t(Throat);
drhodt_average(k)=abs(drhodt(Throat));
dvdt_average(k)=abs(dvdt(Throat));
figure(1)
if k==1
plot(x,mass_flow_rate_1,'b')
hold on
elseif k==50
plot(x,mass_flow_rate_1,'k')
hold on
elseif k==100
plot(x,mass_flow_rate_1,'g')
hold on
elseif k==150
plot(x,mass_flow_rate_1,'y')
hold on
elseif k==200
plot(x,mass_flow_rate_1,'r')
hold on
elseif k==700
plot(x,mass_flow_rate_1,'p')
xlabel('length of the nozzle')
ylabel('mass flow rate')
legend('nt=1','nt=50','nt=100','nt=150','nt=200','nt=700')
grid on
end
end
figure('Name','Timewise variation of density, pressure, Mach no., Temperature at the nozzle throat')
subplot(4,1,1)
plot(rho_th)
xlabel('no. of time step')
ylabel('density')
subplot(4,1,2)
plot(P_th)
xlabel('no. of time step')
ylabel('static pressure')
subplot(4,1,3)
plot(M_th)
xlabel('no. of time step')
ylabel('Mach number')
subplot(4,1,4)
plot(t_th)
xlabel('no. of time step')
ylabel('temperature')
grid on
figure('Name','Timewise variation of the abs values of the time derivatives of rho and V at nozzle throat')
semilogy(drhodt_average)
hold on
semilogy(dvdt_average)
xlabel('number of time step')
ylabel('Residual')
legend('drhodt_{average}','drhodt_{average}')
grid on
end
Calling the function for conservation form:
function [mass_flow_rate_2]=conservative_1D(a,n,x,dx,nt,gamma,Throat)
for i=1:n
if(x(i)>=0 && x(i)<=0.5)
rho(i)=1;
t(i)=1;
elseif(x(i)>=0.5 && x(i)<=1.5)
rho(i)=1.0-0.366*(x(i)-0.5);
t(i)=1.0-0.167*(x(i)-0.5);
elseif(x(i)>=1.5 && x(i)<=3.5)
rho(i)=0.634-0.3879*(x(i)-1.5);
t(i)=0.833-0.3507*(x(i)-1.5);
end
end
v=0.59./(rho.*a);
dt = 0.0267;
U1=rho.*a;
U2=rho.*a.*v;
U3=rho.*(t/(gamma-1)+0.5*(gamma)*v.^2).*a;
for k=1:nt % time loop
rho_old=rho;
v_old=v;
t_old=t;
U1_old=U1;
U2_old=U2;
U3_old=U3;
%pure form of the flux term
F1=U2;
F2=U2.^2./U1+((gamma-1)/gamma).*(U3-0.5*gamma.*(U2.^2./U1));
F3=(gamma.*((U2.*U3)./U1))-(0.5*gamma*(gamma-1).*((U2.^3)./(U1.^2)));
for s=2:n-1
J2(s)=(1/gamma)*rho(s).*t(s).*(a(s+1)-a(s))./dx;
% Predictor method
%continuity equation
dU1dt_p(s)=-(F1(s+1)-F1(s))/dx;
%momentum equation
dU2dt_p(s)=(-(F2(s+1)-F2(s))/dx )+J2(s);
%energy equation
dU3dt_p(s)=-(F3(s+1)-F3(s))/dx;
%updating the values
U1(s)=U1(s)+(dU1dt_p(s))*dt;
U2(s)=U2(s)+(dU2dt_p(s))*dt;
U3(s)=U3(s)+(dU3dt_p(s))*dt;
rho(s)=U1(s)./a(s);
v(s)=U2(s)./U1(s);
t(s)=(gamma-1)*((U3(s)./U1(s))-0.5*gamma*v(s).^2);
end
F1=U2;
F2=U2.^2./U1+((gamma-1)/gamma).*(U3-0.5*gamma.*(U2.^2./U1));
F3=(gamma.*((U2.*U3)./U1))-(0.5*gamma*(gamma-1).*((U2.^3)./(U1.^2)));
for s=2:n-1
J2(s)=(1/gamma)*rho(s).*t(s).*(a(s)-a(s-1))./dx;
% corrector method
dU1dt_c(s)=-(F1(s)-F1(s-1))/dx;
dU2dt_c(s)=(-(F2(s)-F2(s-1))/dx )+J2(s);
dU3dt_c(s)=-(F3(s)-F3(s-1))/dx;
end
%average values
dU1dt=0.5*( dU1dt_p + dU1dt_c );
dU2dt=0.5*( dU2dt_p + dU2dt_c );
dU3dt=0.5*( dU3dt_p + dU3dt_c );
%final value update
for s=2:n-1
U1(s)=U1_old(s)+(dU1dt(s))*dt;
U2(s)=U2_old(s)+(dU2dt(s))*dt;
U3(s)=U3_old(s)+(dU3dt(s))*dt;
end
%applying boundary condition
% AT subsonic iNLET flow
U1(1)=rho(1).*a(1);
U2(1)=2*U2(2)-U2(3);
v(1)=U2(1)/U1(1);
U3(1)=U1(1)*(t(1)/(gamma-1)+gamma*0.5*v(1)^2);
%at supersonic outflow
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);
%decoding the primitive variables
rho=U1./a;
v=U2./U1;
t=(gamma-1)*((U3./U1)-0.5*gamma*v.^2);
%static preesure
P=rho_old.*t_old;
% Mach number
M=v_old./(t_old.^0.5);
% mass flow rate
mass_flow_rate_2=rho_old.*a.*v_old;
% Pressure, temperature, mach number,at the nozzle throat
P_th(k)=P(Throat);
rho_th(k)=rho(Throat);
M_th(k)=M(Throat);
t_th(k)=t(Throat);
dU1dt_avr(k)=abs(dU1dt(Throat));
dU2dt_avr(k)=abs(dU2dt(Throat));
figure(4)
if k==1
plot(x,mass_flow_rate_2,'b')
hold on
elseif k==50
plot(x,mass_flow_rate_2,'k')
hold on
elseif k==100
plot(x,mass_flow_rate_2,'g')
hold on
elseif k==150
plot(x,mass_flow_rate_2,'y')
hold on
elseif k==200
plot(x,mass_flow_rate_2,'r')
hold on
elseif k==700
plot(x,mass_flow_rate_2,'p')
xlabel('length of the nozzle')
ylabel('mass flow rate')
legend('nt=1','nt=50','nt=100','nt=150','nt=200','nt=700')
grid on
end
end
figure('Name','Timewise variation of density, pressure, Mach no. temperature at the nozzle Throat')
subplot(4,1,1)
plot(rho_th)
xlabel('no. of time step')
ylabel('density')
subplot(4,1,2)
plot(P_th)
xlabel('no. of time step')
ylabel('static pressure')
subplot(4,1,3)
plot(M_th)
xlabel('no. of time step')
ylabel('Mach number')
subplot(4,1,4)
plot(t_th)
xlabel('no. of time step')
ylabel('temperature')
grid on
figure('Name','TIMEwise variation of the abs values of the time derivatives of rho and V at nozzle throat')
semilogy(dU1dt_avr)
hold on
semilogy( dU2dt_avr)
xlabel('number of time step')
ylabel('Residual')
legend('dU1dt_{average}',' dU2dt_{average}')
grid on
end
result
Case 1:Non-conservtion form
The above graph shows the variation of p, T, p, and M at the nozzle throat plotted versus the number of time steps. The abscissa starts at zero, which represents the initial conditions, and ends at time step 1000.
We can see largest changes take place at early times, after which the final, steady-state value approached almost asymptotically. Here is the "rubber band effect" at early times the rubber band is "stretched" tightly, and therefore the flow-field variables are driven by a stronger potential and hence change rapidly. At later times, as the steady state is approached, the rubber band is less stretched; it becomes more "relaxed", and the changes become much smaller with time. The horizontal line represent the exact, analytical values as obtained from the equations
The above graph gives the variation of the time derivatives of nondimensional density and velocity as a function of the number of time steps. These are the average time derivatives calculated from Eqs. respectively.
From these results, we can observe two important aspects:
1. At early times, the time derivatives are large, and they oscillate in value. These oscillations are associated with various unsteady compression and expansion waves which propagate through the nozzle during the transient process.
2. At later times, the time derivatives rapidly grow small, changing by six orders of magnitude over a span of 1000 time steps.
The above grapshows the Instantaneous distributions of the nondimensional mass flow as a function of distance through the nozzle at six different times during the time-marching approach to the steady state.
Six different curves are shown, each for a different time during the course of the time-marching procedure. The black coloured curve is the variation of pVA which pertains to the initial conditions(nt=1). The strange-looking, distorted sinelike variation of this curve is simply the product of the assumed initial values for p and V combined with the specified parabolic variation of the nozzle area ratio A. After 50 time steps, the mass flow
distribution through the nozzle has changed considerably; this is given by the curve labeled(nt=50). After 100 time steps , the mass flow distribution has changed radically; the mass flow variation is simply flopping around inside the nozzle due to the transient variation of the flow-field variables. However, after 200 time steps, the mass flow distribution is beginning to settle down, and after 700 time
steps , the mass flow distribution is a straight, horizontal line across the graph. This says that the mass flow has converged to a constant, steady-state value throughout the nozzle.
Case 2: Conservation form
The above graph shows the variation of p, T, p, and M at the nozzle throat plotted versus the number of time steps for conservation form
The above graph shows the variation of abs values of the time derivative of solution vectors U1 and U2
At early times, the time derivatives are large, and they oscillate in value. These oscillations are associated with various unsteady compression and expansion waves which propagate through the nozzle during the transient process. At later times, the time derivatives rapidly grow small, changing by six orders of magnitude over a span of 1000 time steps.
The above graph shows Variation of mass flow distribution through the nozzle at different times during the time-marching process, solution of the conservation form of the governing equations. Star curve indicates the exact analytical value.
Variation of mass flow distribution through the nozzle at different times for both non conversation and conservation form is almost same but with aslight difference .
The above graph shows the detailed comparison of the steady-state mass flow variations obtained with the nonconservation and conservation forms of the governing flow equations.
Here we see that the steady-state mass flow distribution with the predicted by the conservation form of the equations is much more satisfactory than that obtained from the nonconservation form, on two accounts:
- The conservation form gives a distribution that is much closer to being a constant. In contrast, the nonconservation results have (on the magnified scale) a sizeable _variation, with some spurious oscillations at both the inflow and outflow boundaries.
- The steady-state mass flow results obtained with the conservation form are, on the whole, much closer to the exact analytical solution of p'A'V' = 0.579, shown by star line.
Grid Independence test:
| `rho^*/rho^0` | `t^*/t^0` | `p^*/p^0` | M | |
| Exact analytical solution | 0.634 | 0.833 | 0.528 | 1.00 |
| NOn conservation form(n=31) | 0.6388 | 0.8365 | 0.5343 | 0.9993 |
| Conservation form(n=31) | 0.6486 | 0.834 | 0.5281 | 1.0196 |
| Coservation form (n=61) | 0.6504 | 0.8415 | 0.5281 | 1.0196 |
Conclusion:
- The conservation form yields a better mass flow distribution. The conservation form simply does a better job of conserving mass.
- The nonconservation form leads to smaller residuals. The amount by which the residuals decay is often used as an index of "quality" of the numerical algorithm. In this sense, the nonconservation form does a better job.
- 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.
- 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.