Courses by Software
Courses by Semester
Courses by Domain
Tool-focused Courses
Machine learning
POPULAR COURSES
Post Graduate Program in Hybrid Electric Vehicle Design and AnalysisPost Graduate Program in Computational Fluid DynamicsPost Graduate Program in CADPost Graduate Program in CAEPost Graduate Program in Manufacturing DesignPost Graduate Program in Computational Design and Pre-processingPost Graduate Program in Complete Passenger Car Design & Product DevelopmentSuccess Stories
NUMERICAL SOLUTIONS OF QUASI-ONE-DIMENSIONAL NOZZLE FLOWS USING MACCORMACK METHOD
Aim: To simulate the isentropic flow through a quasi 1D subsonic-supersonic nozzle using MacCormack's technique. Abstract: The flow is steady, isentropic, through a convergent-divergent nozzle as shown in the figure below. The flow at the inlet of the nozzle comes from the reservoir where the pressure and temperature…
Faizan Akhtar
updated on 04 Apr 2021
Aim: To simulate the isentropic flow through a quasi 1D subsonic-supersonic nozzle using MacCormack's technique.
Abstract: The flow is steady, isentropic, through a convergent-divergent nozzle as shown in the figure below. The flow at the inlet of the nozzle comes from the reservoir where the pressure and temperature are denoted by p0p0 and T0T0 respectively. The cross-sectional of the reservoir is very large (theoretically A→⋈A→⋈) and hence the velocity is very small (V==0V==0). The flow expands isentropically at the supersonic speed at the exit, where the exit pressure, temperature, velocity, and Mach number are denoted by pepe, TeTe, VeVe, MeMe. The flow is subsonic (M<1M<1) at the inlet of the nozzle, sonic at the throat section (M=1M=1).
Please note that in the energy equation e has been replaced with T
Non-conservative equation
The integral form of the continuity equation is represented below
∂∂t∫∫V∫ρ∗dV+∫∫Sρ∗v∗dS=0∂∂t∫∫V∫ρ∗dV+∫∫Sρ∗v∗dS=0
The equation is applied to the shaded control volume
The left side of the figure consist of pp, vv,ρρ,TT and the right side of the variable is represented as p+dpp+dp, v+dvv+dv, ρ+dρρ+dρ,T+dTT+dT. They are uniform over the area A+dAA+dA.
Continuity equation
The volume integral is represented as
∂∂t∫∫V∫ρ∗dV=∂∂t(ρ∗A∗dx)∂∂t∫∫V∫ρ∗dV=∂∂t(ρ∗A∗dx)
The surface integral is represented as
∫∫Sρ∗v∗dS=-ρ∗v∗A+(ρ+dρ)(v+dv)(A+dA)∫∫Sρ∗v∗dS=−ρ∗v∗A+(ρ+dρ)(v+dv)(A+dA)
The -− sign in the equation is because dS always point outside the control volume.
Thus expanding the above equation
∫∫Sρ∗v∗dS=-ρ∗v∗A+ρ∗v∗A+ρ∗v∗dA∫∫Sρ∗v∗dS=−ρ∗v∗A+ρ∗v∗A+ρ∗v∗dA
+ρ∗A∗dv+ρ∗dv∗dA+A∗v∗dρ+v∗dA∗dρ+A∗dv∗dρ+dρ∗dv∗dA+ρ∗A∗dv+ρ∗dv∗dA+A∗v∗dρ+v∗dA∗dρ+A∗dv∗dρ+dρ∗dv∗dA
The term consisting of product of differential approaches to zero. Thus the equation becomes
∫∫Sρ∗v∗dS=ρ∗v∗dA+ρ∗A∗dv+A∗v∗dρ=d(ρ∗A∗v)∫∫Sρ∗v∗dS=ρ∗v∗dA+ρ∗A∗dv+A∗v∗dρ=d(ρ∗A∗v)
Thus the equation becomes
∂∂t(ρ∗A∗dx)+d(ρ∗A∗v)=0∂∂t(ρ∗A∗dx)+d(ρ∗A∗v)=0
Thus dividing by dx throughout
∂∂t(ρ∗A)+∂∂x(ρ∗A∗v)=0∂∂t(ρ∗A)+∂∂x(ρ∗A∗v)=0
This equation represent the continuity equation for unsteady,quasi one-dimensional flow.
Momentum equation
Writing the integral form of the equation as
∂∂t∫∫V∫(ρ∗u)∗dV+∫∫s(ρ∗u∗V)∗dS=-∫∫s(p∗dS)x∂∂t∫∫V∫(ρ∗u)∗dV+∫∫s(ρ∗u∗V)∗dS=−∫∫s(p∗dS)x
∂∂t∫∫V∫(ρ∗u)∗dV=∂∂t(ρ∗v∗A∗dx)∂∂t∫∫V∫(ρ∗u)∗dV=∂∂t(ρ∗v∗A∗dx)
∫∫s(ρ∗u∗V)∗dS=-ρ∗v2∗A+(ρ+dρ)(v+dv)2(A+dA)∫∫s(ρ∗u∗V)∗dS=−ρ∗v2∗A+(ρ+dρ)(v+dv)2(A+dA)
∫∫s(p∗dS)x=-p∗A+(p+dp)∗(A+dA)-2∗p∗(dA2)∫∫s(p∗dS)x=−p∗A+(p+dp)∗(A+dA)−2∗p∗(dA2)
Solving for each term and cancelling the product of differential terms yield
∂∂t(ρ∗v∗A∗dx)+d(ρ∗v2∗A)=-Adp∂∂t(ρ∗v∗A∗dx)+d(ρ∗v2∗A)=−Adp
Dividing by dx and taking limit dx→0dx→0 yields partial differential equation of the form
∂∂t(ρ∗v∗A)+∂∂x(ρ∗v2∗A)=-A∗∂p∂x∂∂t(ρ∗v∗A)+∂∂x(ρ∗v2∗A)=−A∗∂p∂x
This equation represent conservation form of the momentum equation.
Taking the vv term out of the differential equation gives
∂∂t(ρ∗v∗A)+∂∂x(ρ∗v2∗A)=-A∗∂p∂x∂∂t(ρ∗v∗A)+∂∂x(ρ∗v2∗A)=−A∗∂p∂x
Multiplying continuity equation by v
v∗∂∂t(ρ∗A)+v∗∂∂x(ρ∗A∗v)=0v∗∂∂t(ρ∗A)+v∗∂∂x(ρ∗A∗v)=0
Substracting the equation
∂∂t(ρ∗v∗A)-v∗∂∂t(ρ∗A)+∂∂x(ρ∗v2∗A)-v∗∂∂x(ρ∗A∗v)=-A∗∂p∂x∂∂t(ρ∗v∗A)−v∗∂∂t(ρ∗A)+∂∂x(ρ∗v2∗A)−v∗∂∂x(ρ∗A∗v)=−A∗∂p∂x
Expanding the derivatives and canceling the like terms
ρ∗A∗∂v∂t+ρ∗A∗v∗∂v∂x=-A∗∂p∂xρ∗A∗∂v∂t+ρ∗A∗v∗∂v∂x=−A∗∂p∂x
Dividing by A
ρ∗∂v∂t+ρ∗v∗∂v∂x=-∂p∂xρ∗∂v∂t+ρ∗v∗∂v∂x=−∂p∂x
The equation represent momentum equation for quasi-1d flow in non conservative form.
Energy equation
For an adiabatic flow (q=0q=0) with no body forces and no viscous effects, the integral form of the equation is represented as
∂∂t∫∫V∫(e+v22)dV+∫∫S(ρ+v22)∗V∗dS=-∫∫s(ρ∗V)dS∂∂t∫∫V∫(e+v22)dV+∫∫S(ρ+v22)∗V∗dS=−∫∫s(ρ∗V)dS
Expanding the above equation
∂∂t[ρ∗(e+v22)∗A∗dx]-ρ∗(e+v22)∗V∗A+(ρ+dρ)∗[e+de+(v+dv)22]∗(v+dv)∗(A+dA)∂∂t[ρ∗(e+v22)∗A∗dx]−ρ∗(e+v22)∗V∗A+(ρ+dρ)∗[e+de+(v+dv)22]∗(v+dv)∗(A+dA)
=-[-ρ∗v∗A+(p+dp)∗(v+dv)∗(A+dA)-2∗(p∗v∗dA2)]=−[−ρ∗v∗A+(p+dp)∗(v+dv)∗(A+dA)−2∗(p∗v∗dA2)]
Neglecting the product of differentiation and cancelling the like terms we have
∂∂t[ρ(e+v22)∗A∗dx]+d(ρ∗e∗v∗A)+d∗(ρ∗v3∗A)2=-d(ρ∗A∗v)∂∂t[ρ(e+v22)∗A∗dx]+d(ρ∗e∗v∗A)+d∗(ρ∗v3∗A)2=−d(ρ∗A∗v)
or
∂∂t[ρ(e+v22)∗A∗dx]+d[ρ(e+v22)∗v∗A]=-d(ρ∗A∗v)∂∂t[ρ(e+v22)∗A∗dx]+d[ρ(e+v22)∗v∗A]=−d(ρ∗A∗v)
Dividing by dx throughout and taking limit dx→0dx→0 the partial differential equation is obtained
∂(ρ∗(e+v22)∗A)∂t+∂(ρ∗(e+v22))∂x⋅=-∂(ρ∗A∗v)∂x∂(ρ∗(e+v22)∗A)∂t+∂(ρ∗(e+v22))∂x⋅=−∂(ρ∗A∗v)∂x
The above equation is the conservation form of the energy equation.
In order to achieve non conservative form of energy equation
Multiplying v in conservation form of the momentum equation
∂∂t(ρ∗v22∗A)+∂∂x(ρ∗v33∗A)=-A∗v∗∂p∂x∂∂t(ρ∗v22∗A)+∂∂x(ρ∗v33∗A)=−A∗v∗∂p∂x
Substracting from conservation form of energy equation
∂(ρ∗e∗A)∂t+∂(ρ∗e∗v∗A)∂x=-p∗∂(A∗v)∂x∂(ρ∗e∗A)∂t+∂(ρ∗e∗v∗A)∂x=−p∗∂(A∗v)∂x
Multiplying the above equation by e gives
e∗∂(ρ∗A)∂t+e∗∂(ρ∗A∗v)∂x=0e∗∂(ρ∗A)∂t+e∗∂(ρ∗A∗v)∂x=0
ρ∗A∗∂e∂t+(ρ∗A∗v)∗∂e∂x=-p∂(A∗v)∂xρ∗A∗∂e∂t+(ρ∗A∗v)∗∂e∂x=−p∂(A∗v)∂x
Dividing by A throughout
ρ∗∂e∂t+ρ∗v∗∂e∂x=-ρ∗∂v∂x-ρ∗v∗∂(lnA)∂xρ∗∂e∂t+ρ∗v∗∂e∂x=−ρ∗∂v∂x−ρ∗v∗∂(lnA)∂x
The above equation is non-conservative form of the energy equation.
Non dimensional form of the equation
Please note that while using the Maccormack method we will be using the non-dimensional form of the above equation for ease in the calculation.
Non-dimensional temperature T′=TT0T'=TT0 T0T0 is the reservoir temperature
Non-dimensional density ρ′=ρρ0ρ'=ρρ0 ρ0ρ0 is the reservoir density.
Non-dimensional length x′=xLx'=xL LL is the length of nozzle.
Non-dimensional velocity v′=va0v'=va0 a0a0 is the speed of the sound in the reservoir.
Non-dimensional time t′=tLa0t'=tLa0
Non-dimensional area A′=AA∗A'=AA∗ A∗A∗ is the sonic throat area.
Replacing the above equation with non-dimensional variables we get
Continuity equation = ∂ρ′∂t′=-ρ′∗∂v′∂x′-ρ′∗v′∗∂(lnA′)∂x′-v′∗∂ρ′∂x′
Momentum equation= ∂v′∂t′=-v′∗∂v′∂x′-1γ∗(∂T′∂x′+T′ρ′∗∂ρ′∂x′)
Energy equation ∂T′∂t′=-v′∗∂T′∂x′-(γ-1)∗T′∗[∂v′∂x′+v′∗∂(lnA′)∂x′]
Conservative equation
The conservative form of equation can be expressed in terms of solution vector U flux vector F and the source term J which are as follows
U1=ρ′∗A′
U2=ρ′∗A′∗v′
U3=ρ′(T′γ-1+γ2∗v′2)∗A′
F1=ρ′∗A′∗v′
F2=U22U1+γ-1γ(U3U1-γ2∗(U2U1)2)
F3=γ∗(U2∗U3U1)-γ(γ-1)2∗U32U21
J2=ρ′∗T′γ∗∂A∂x
MacCormack method
It is used for discretising hyperbolic partial differential equation. It involves predictor and corrector method.
Predictor method: It uses forward differencing for space derivative. In this the provisional values of variables (ρ Tv) at (n+1) timesteps.
Corrector method: It uses backward differencing for space derivative. In this the predicted values are corrected. The timesteps are used △t2
The method is used for non linear equation (inviscid Burgers equation, Euler equation). For non linear equation it provides the best result. For linear equation Lax-Wendroff method is used.
The above concept have been applied to produce the Matlab code as shown below
- Matlab code for the main program.
% main program
clear
close all
clc
% defining one dimensional mesh size
L=3;% length of nozzle
n=31;% number of grid points
x=linspace(0,L,n);
dx=x(2)-x(1);
% defining other input parameters
gamma=1.4;
C=0.5;
nt=linspace(1,1400,1400);
% writing function for non conservative form
[mfr_analytical,mfr_e]=nozzle(x,L,n,gamma,nt);
[mfr_con]=conservative(x,L,n,gamma,nt);
figure(7);
plot(x,mfr_analytical)
hold on
plot(x,mfr_con,'k')
hold on
plot(x,mfr_e,'r')
xlabel('Length of nozzle')
ylabel('Non dimensional mass flow rate')
legend('Analytical mass flow rate','Mass flow rate for conservative form','Mass flow rate for non conservative form')
title('Comparison of Normalised mass flow rate distributions of both the forms')
- Creating user-defined function for the conservative equation.
function [mfr_con]=conservative(x,L,n,gamma,nt)
%% Conservative form
clear
close all
clc
% defining one dimensional mesh size
L=3;% length of nozzle
n=31;% number of grid points
x=linspace(0,L,n);
dx=x(2)-x(1);
a=1+2.2*(x-1.5).^2;% area
% throat
throat=find(a==1);
%CFL number
C=0.5;
% defining other input parameters
gamma=1.4;
% intialising temperature and density
for i=1:length(x)
if x(i)>=0 && x(i)<0.5
rho_con(i)=1;
temp_con(i)=1;
elseif x(i)>=0.5 && x(i)<1.5
rho_con(i)=1-(0.366.*(x(i)-0.5));
temp_con(i)=1-(0.167.*(x(i)-0.5));
elseif x(i)>=1.5 && x(i)<=3
rho_con(i)=0.634-(0.3879.*(x(i)-1.5));
temp_con(i)=0.833-(0.3507.*(x(i)-1.5));
end
end
% calculating primitive variables
v_con=0.59./(rho_con.*a);
p_con=rho_con.*temp_con;
% calculating solution vector
U1=rho_con.*a;
U2=rho_con.*a.*v_con;
U3=rho_con.*a.*((temp_con)/(gamma-1)+(0.5*gamma*v_con.^2));
nt=linspace(1,1400,1400);
tic;
iteration=0;
for k=1:length(nt)
% predictor step
% creating copy of solution vector
U1_old=U1;
U2_old=U2;
U3_old=U3;
% calculating dt
dt_con=min((C*dx)./((sqrt(temp_con)+v_con)));
% calculating flux vector
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);
% calculating source term
for j=2:n-1
J2(j)=(rho_con(j).*temp_con(j).*(a(j+1)-a(j)))/(gamma*dx);
dU1dt_p(j)=-(F1(j+1)-F1(j))/dx;
dU2dt_p(j)=-(F2(j+1)-F2(j))/(dx)+J2(j);
dU3dt_p(j)=-(F3(j+1)-F3(j))/dx;
% updating values
U1(j)=U1(j)+(dU1dt_p(j)*dt_con);
U2(j)=U2(j)+(dU2dt_p(j)*dt_con);
U3(j)=U3(j)+(dU3dt_p(j)*dt_con);
end
% updating primitive values
rho_con=U1./a;
v_con=U2./U1;
temp_con=(gamma-1)*(U3./U1-(gamma/2)*v_con.^2);
p_con=rho_con.*temp_con;
% updating flux vector
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
for j=2:n-1
J2(j)=(rho_con(j).*temp_con(j).*(a(j)-a(j-1)))/(gamma*dx);
dU1dt_c(j)=-(F1(j)-F1(j-1))/dx;
dU2dt_c(j)=-(F2(j)-F2(j-1))/(dx)+J2(j);
dU3dt_c(j)=-(F3(j)-F3(j-1))/dx;
end
% average time step value
dU1dt_a=0.5*(dU1dt_p+dU1dt_c);
dU2dt_a=0.5*(dU2dt_p+dU2dt_c);
dU3dt_a=0.5*(dU3dt_p+dU3dt_c);
for j=2:n-1
U1(j)=U1_old(j)+(dU1dt_a(j)*dt_con);
U2(j)=U2_old(j)+(dU2dt_a(j)*dt_con);
U3(j)=U3_old(j)+(dU3dt_a(j)*dt_con);
end
iteration=iteration+1;
% left boundary condition
U1(1)=rho_con(1).*a(1);
U2(1)=2*U2(2)-U2(3);
U3(1)=U1(1)*(temp_con(1)./(gamma-1)+(gamma/2)*v_con(1).^2);
% right boundary condition
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);
% calculating the primitive variables at exit
rho_con=U1./a;
v_con=U2./U1;
temp_con=(gamma-1)*(U3./U1-(gamma/2)*v_con.^2);
p_con=rho_con.*temp_con;
mfr_con=rho_con.*v_con.*a;
M_con=v_con./temp_con.^0.5;
% calculating throat values
rho_con_throat(k)=rho_con(throat);
v_con_throat(k)=v_con(throat);
temp_con_throat(k)=temp_con(throat);
mfr_con_throat(k)=mfr_con(throat);
M_con_throat(k)=M_con(throat);
p_throat(k)=p_con(throat);
mfr(iteration,:)=rho_con.*a.*v_con;
end
computational_time=toc;
figure(4);
plot(x,mfr(1,:))
hold on
plot(x,mfr(50,:))
hold on
plot(x,mfr(100,:))
hold on
plot(x,mfr(150,:))
hold on
plot(x,mfr(200,:))
hold on
plot(x,mfr(700,:))
hold on
plot(x,mfr(1400,:))
hold on
xlabel('Length of nozzle')
ylabel('Non dimensional mass flow rate')
legend('iteration=1','iteration=50','iteration=100''iteration=150','iteration=200','iteration=700','iteration=1400')
title('Variation of non dimensional mass flow rate in conservative form with time steps=1400.')
figure(5);
subplot(5,1,1)
plot(nt,rho_con_throat)
xlabel('Number of time steps')
ylabel('Density')
legend('Density')
title('Timewise variation of primitive variables at throat for timesteps 1400 in conservative form')
subplot(5,1,2)
plot(nt,v_con_throat)
xlabel('Number of time steps')
ylabel('Velocity')
legend('Velocity')
subplot(5,1,3)
plot(nt,temp_con_throat)
xlabel('Number of time steps')
ylabel('Temperature')
legend('Temperature')
subplot(5,1,4)
plot(nt,M_con_throat)
xlabel('Number of time steps')
ylabel('Mach number')
legend('Mach number')
subplot(5,1,5)
plot(nt,p_throat)
xlabel('Number of time steps')
ylabel('Pressure')
legend('Pressure')
figure(6);
subplot(5,1,1)
plot(x,rho_con)
xlabel('Length of nozzle')
ylabel('Density')
legend('Density')
title('Length wise variation of primitive variables for conservative form')
subplot(5,1,2)
plot(x,v_con)
xlabel('Length of nozzle')
ylabel('Velocity')
legend('Velocity')
subplot(5,1,3)
plot(x,temp_con)
xlabel('Length of nozzle')
ylabel('Temperature')
legend('Temperature')
subplot(5,1,4)
plot(x,p_con)
xlabel('Length of nozzle')
ylabel('Pressure')
legend('Pressure')
subplot(5,1,5)
plot(x,M_con)
xlabel('Length of nozzle')
ylabel('Mach number')
legend('Mach number')
end- Creating user-defined function for the non-conservative equation.
function [mfr_analytical,mfr_e]=nozzle(x,L,n,gamma,nt)
%Non conservative form
clear
close all
clc
% defining one dimensional mesh size
L=3;% length of nozzle
n=31;% number of grid points
x=linspace(0,L,n);
dx=x(2)-x(1);
% defining other input parameters
gamma=1.4;
C=0.5;
nt=linspace(1,1400,1400);
% nozzle properties at initial condition
mfr_analytical=0.579*ones(1,n);
a=1+2.2*(x-1.5).^2;% area
rho=1-0.3146.*x;% density
t=1-0.2314.*x;% temperature
v=(0.1+1.09*x).*t.^0.5;% velocity
mfr_initial=rho.*a.*v; % mass flow rate initial
M_initial=v./t.^0.5; % mach number initial
p_initial=rho.*t; % pressure initial
% calculating throat value
throat=find(a==1);
% time interval value
dt=min((C*dx)./((sqrt(t)+v)));
number_iter=0;
% starting time loop
tic;
for k=1:length(nt)
% creating a copy of profile for comparison
rho_old=rho;
v_old=v;
t_old=t;
% space loop
for j=2:n-1
% defining terms for calculation convenience
% predictor step
dvdx_p=(v(j+1)-v(j))/dx;
dlogadx_p=(log(a(j+1))-log(a(j)))/dx;
drhodx_p=(rho(j+1)-rho(j))/dx;
dtdx_p=(t(j+1)-t(j))/dx;
% continuity equation
drhodt_p(j)=(-rho(j).*dvdx_p)+(-rho(j).*v(j).*dlogadx_p)+(-v(j).*drhodx_p);
% momentum equation
dvdt_p(j)=(-v(j).*dvdx_p)+((-1/gamma)*(dtdx_p+(t(j)./rho(j)).*(drhodx_p)));
% energy equation
dtdt_p(j)=(-v(j).*dtdx_p)-(gamma-1).*t(j).*(dvdx_p+(v(j).*dlogadx_p));
% updating v,rho and t
v(j)=v(j)+dvdt_p(j).*dt;
rho(j)=rho(j)+drhodt_p(j).*dt;
t(j)=t(j)+dtdt_p(j).*dt;
end
for j=2:n-1
% corrector step
dvdx_c=(v(j)-v(j-1))/dx;
dlogadx_c=(log(a(j))-log(a(j-1)))/dx;
drhodx_c=(rho(j)-rho(j-1))/dx;
dtdx_c=(t(j)-t(j-1))/dx;
% continuity equation
drhodt_c(j)=(-rho(j).*dvdx_c)+(-rho(j).*v(j).*dlogadx_c)+(-v(j).*drhodx_c);
% momentum equation
dvdt_c(j)=(-v(j).*dvdx_c)+((-1/gamma)*(dtdx_c+(t(j)./rho(j)).*(drhodx_c)));
% energy equation
dtdt_c(j)=(-v(j).*dtdx_c)-(gamma-1).*t(j).*(dvdx_c+(v(j).*dlogadx_c));
end
% average time step
dvdt_a=0.5*(dvdt_c+dvdt_p);
dtdt_a=0.5*(dtdt_c+dtdt_p);
drhodt_a=0.5*(drhodt_c+drhodt_p);
% updating final values
for j=2:n-1
v(j)=v_old(j)+dvdt_a(j).*dt;
t(j)=t_old(j)+dtdt_a(j).*dt;
rho(j)=rho_old(j)+drhodt_a(j).*dt;
end
number_iter=number_iter+1;
% applying boundary condition
v(1)=2*v(2)-v(3);
rho(1)=1;
t(1)=1;
v(n)=2*v(n-1)-v(n-2);
t(n)=2*t(n-1)-t(n-2);
rho(n)=2*rho(n-1)-rho(n-2);
% primitive variables at exit
p_e=rho.*t;
M_e=v./t.^0.5;
mfr_e=rho.*a.*v;
% value at throat values
mfr_throat(k)=mfr_e(throat);
p_throat(k)=p_e(throat);
M_throat(k)=M_e(throat);
v_throat(k)=v(throat);
rho_throat(k)=rho(throat);
t_throat(k)=t(throat);
figure(1);
if (number_iter==1)
plot(x,mfr_e)
hold on
elseif (number_iter==50)
plot(x,mfr_e)
hold on
elseif (number_iter==100)
plot(x,mfr_e)
hold on
elseif (number_iter==150)
plot(x,mfr_e)
hold on
elseif (number_iter==200)
plot(x,mfr_e)
hold on
elseif (number_iter==700)
plot(x,mfr_e)
hold on
elseif (number_iter==1400)
plot(x,mfr_e)
hold on
xlabel('Length of nozzle')
ylabel('Non dimensional mass flow rate')
legend('number iter==1','number iter==50','number iter==100','number iter==200','number iter==700','number iter==1400','mfr analytical','mfr initial')
title('Variation of non dimensional mass flow rate along the length of nozzle with varying time steps in NC form.')
end
end
computational_time=toc;
figure(2);
subplot(5,1,1)
plot(nt,p_throat)
xlabel('Number of time steps')
ylabel('Pressure')
legend('Pressure at throat')
title('Variation of primitive variables at throat for time steps 1400 in NC form')
subplot(5,1,2)
plot(nt,rho_throat)
xlabel('Number of time steps')
ylabel('Density')
legend('Density at throat')
subplot(5,1,3)
plot(nt,t_throat)
xlabel('Number of time steps')
ylabel('Temperature')
legend('Temperature at throat')
subplot(5,1,4)
plot(nt,v_throat)
xlabel('Number of time steps')
ylabel('Velocity')
legend('Velocity at throat')
subplot(5,1,5)
plot(nt,M_throat)
xlabel('Number of time steps')
ylabel('Mach number')
legend('Mach number at throat')
figure(3);
subplot(5,1,1)
plot(x,p_e)
xlabel('Length of nozzle')
ylabel('Pressure')
legend('Pressure at exit')
title('Variation of primitive variables at exit along the length of nozzle')
subplot(5,1,2)
plot(x,rho)
xlabel('Length of nozzle')
ylabel('Density')
legend('Density at exit')
subplot(5,1,3)
plot(x,t)
xlabel('Length of nozzle')
ylabel('Temperature')
legend('Temperature at exit')
subplot(5,1,4)
plot(x,v)
xlabel('Length of nozzle')
ylabel('Velocity')
legend('Velocity at exit')
subplot(5,1,5)
plot(x,M_e)
xlabel('Length of nozzle')
ylabel('Mach number')
legend('Mach number at exit')
endDemonstration of grid independence.
| S No |
Number of grid points. |
Condition at the nozzle throat | |||||||
| Conservative form | Non Conservative form | ||||||||
|
Density ρ |
Temperature T |
Pressure p |
Mach number M |
Density ρ |
Temperature T |
Pressure p |
Mach number M |
||
| 1 | 31 | 0.6504 | 0.8400 | 0.5464 | 0.9817 | 0.6387 | 0.8365 | 0.5342 | 0.9994 |
| 2 | 93 | 0.6373 | 0.8359 | 0.5327 | 1.0001 | 0.6399 | 0.8365 | 0.5353 | 1.0035 |
| 3 | 155 | 0.6415 | 0.8372 | 0.5371 | 1.0034 | 0.6666 | 0.8507 | 0.5671 | 1.0197 |
| Exact analytical solution | 0.634 | 0.833 | 0.528 | 1.000 | 0.634 | 0.833 | 0.528 | 1.000 | |
The above table represents the steady state solution at the throat section of nozzle for both the cases. The same is true for all the location within the nozzle. The deviation from the exact solution is very marginal for grid number 93 and 155 respectively. Thus increasing the grid points will increase the numerical accuracy of the solution, but it will also increase the execution time which is to be kept in mind.
Thus we can say that instead of increasing the grid points we can execute the simulation for 31 grid points which will save our time.
Courant number effect
| S No | Courant Number | Condition at the nozzle throat | |||||||
| Conservative form | Non Conservative form | ||||||||
|
Density ρ |
Temperature T |
Pressure p |
Mach number M |
Density ρ |
Temperature T |
Pressure p |
Mach number M |
||
| 1 | 0.5 | 0.6504 | 0.8400 | 0.5464 | 0.9817 | 0.6387 | 0.8365 | 0.5342 | 0.9994 |
| 2 | 0.6 | 0.6478 | 0.8387 | 0.5433 | 0.9867 | 0.6389 | 0.8365 | 0.5344 | 0.9993 |
| 3 | 0.7 | 0.6457 | 0.8376 | 0.5408 | 0.9906 | 0.6391 | 0.8366 | 0.5346 | 0.9992 |
| 4 | 0.8 | 0.6441 | 0.8368 | 0.5389 | 0.9936 | 0.6393 | 0.8367 | 0.5349 | 0.9991 |
| 5 | 0.9 | 0.6430 | 0.8362 | 0.5376 | 0.9957 | 0.6395 | 0.8368 | 0.5351 | 0.9989 |
| 6 | 1.2 | Solution did not converge/Program is unstable and it blew up. | |||||||
| Exact analytical sloution | 0.634 | 0.833 | 0.528 | 1.000 | 0.634 | 0.833 | 0.528 | 1.000 | |
From the above table it can be concluded that when courant number is increased to 1.2 instabilities do occur and the program blows up. The CFL number holds good for flow involving linear equations, but for non-linear parabolic partial differential equation it can be useful for providing the values of △t.
Expected graphs
- Steady-state distribution of primitive variables inside the nozzle for the conservative equation.
- Steady-state distribution of primitive variables inside the nozzle for the non-conservative equation.
- Time-wise variation of the primitive variables for the conservative equation.
- Time-wise variation of the primitive variables for the non-conservative equation.
- Variation of Mass flow rate distribution inside the nozzle at different time steps during the time-marching process for the conservative equation.
- Variation of Mass flow rate distribution inside the nozzle at different time steps during the time-marching process for the non-conservative equation.
- Comparison of Normalized mass flow rate distributions of both forms
Conclusion
- The program was written for both the forms and was executed in matlab.
- It can be inferred from the "steady state distribution of primitive variables inside the nozzle for both the forms'' that there is no much difference of the values in both the form as it can be observed in grid independence test.
- It can be inferred from the "time wise variation of the primitive variables for both the forms" that the largest changes takes place at the earliest time-steps, after 400 timesteps the value of the primitive variables reaches to the steady state.
- It can be inferred from the "variation of the non dimensional mass flow rate along the nozzle for different time steps in conservative form" that the result after 100 iteration shows a humped distribution. The mass flow rate has started to converge after 200 iteration and after 700 iteration it has started to converge and becomes equal to 0.579.
- It can be inferred from the "variation of the non dimensional mass flow rate along the nozzle for different time steps in non conservative form" that the first iteration value is represented by the blue line graph. After 50 iteration the mass flow rate has changed drastically, and after 100 iteration the mass flow rate has changed radically. After 700 iteration the mass flow rate begins to settle down, and in between 700 to 1400 iteration the mass flow rate has converged and becomes steady.
- It can be inferred from the " comparison of normalised mass flow rate distribution of both the forms" that the mass flow rate predicted by the conservation form of the equation is much more satisfactory than the one obtained from the nonconservation form of the equation on the basis of two accounts: that the conservation form is much closer to 0.579, on the other hand in non-conservative equations some spurious oscillations can be found at the inlet and outlet of the nozzle. The conservation form of the equation is closer to analytical mass flow rate.
- The execution time for the conservation form of equation is 1.3661 and for the non-conservative form is 2.4725.
Source: CFD book by John Anderson, Google
Leave a comment
Thanks for choosing to leave a comment. Please keep in mind that all the comments are moderated as per our comment policy, and your email will not be published for privacy reasons. Please leave a personal & meaningful conversation.
Other comments...
Be the first to add a comment
Read more Projects by Faizan Akhtar (34)
Week-4 : Basic Calibration of Single cylinder CI-Engine
Aim: Basic Calibration of Single cylinder CI-Engine Objective : Explore Tutorial No 1- CI final 1.Compare SI vs CI and list down differences (assignment no 2-SI) 2. Comments on the following parameters BSFC Exhaust Temperature A/F ratios 3.Change MFB 50 and observe impact on performance Introduction Difference…
11 Nov 2021 05:26 AM IST
Week 2 : Basic Calibration of Single cylinder SI-Engine
Aim: Basic Calibration of Single-cylinder SI-Engine Objective: 1. Run the case at 1800 rpm and list down important parameters (20 Marks) air flow rate BMEP BSFC In-cylinder pressure 2. Increase the power output at 3600 rpm by 10% (30 Marks) Introduction A spark-ignition engine (SI engine) is…
22 Oct 2021 07:11 AM IST
Week 1 : Exploring the GUI of GT-POWER
Aim: Exploring the GUI of GT-suite GT-suite can be used in a wide range of applications. There are several industries that are using GT-suite for various applications On-highway and off-highway vehicle Marine and rail Industrial machinery Aerospace Power generation Multiphysics platform GT-Suite has a versatile multiphysics…
12 Oct 2021 02:52 PM IST
Week 8: Literature review - RANS derivation and analysis
Aim: Literature review - RANS derivation and analysis Objective: Apply Reynolds decomposition to the NS equations and come up with the expression for Reynolds stress. Explain your understanding of the terms Reynolds stress What is turbulent viscosity? How is it different from molecular viscosity? Introduction…
01 Sep 2021 07:52 AM IST
Related Courses
Skill-Lync offers industry relevant advanced engineering courses for engineering students by partnering with industry experts.
Our Company
4th Floor, BLOCK-B, Velachery - Tambaram Main Rd, Ram Nagar South, Madipakkam, Chennai, Tamil Nadu 600042.
Top Individual Courses
Top PG Programs
Skill-Lync Plus
Trending Blogs
© 2025 Skill-Lync Inc. All Rights Reserved.