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
1D Super-sonic nozzle flow simulation using Macormack Method in Octave / MATLAB
Aim:-> To perform numerical flow simulation of quasi 1D subsonic-supersonic nozzle using MacCormack technique with conservative and non conservative methods. Theory:- MacCormack technique:- MacCormack method is an explicit finite difference technique which is second order accurate in both time…
Smitesh Badnapurkar
updated on 26 Jan 2021
Aim:-> To perform numerical flow simulation of quasi 1D subsonic-supersonic nozzle using MacCormack technique with conservative and non conservative methods.
Theory:-
MacCormack technique:- MacCormack method is an explicit finite difference technique which is second order accurate in both time and space. The method is applied for discretizing hyperbolic partial differential equation. As in MacCormack method there are two steps i.e. predictor step with forward differencing and corrector step with backward differencing, the solution obtained by taking average of predictor and corrector step gives second order accurate. The procedure for this technique is as follows.
- calculate initial values for equation
- To calculate predictor step values by using forward differencing for space terms.
- Update values on RH side of equation
- To calculate corrector step values by using backward differencing for space term.
- Take average of corrector step values and predictor step values. Corrector value+Predictor value2Corrector value+Predictor value2
- Update RH side values by taking previous timestep values
- Repeat above all steps for successive timestep until solution gets converged.
There are two methods to solve almost every CFD problem numerically i.e. by using conservative equations and non conservative equations. From theoretical and mathematical view both forms are equivalent and can be converted into one another. For the engineering problems, user can choose any of the above flow method, each method possess its own advantages and disadvantages.
~ Conservative form
Consider an infinitesimal small fluid element inside nozzle, and the fluid moving through that element, such type of condition is conservative form and the equation obtained in such cases are called as conservative form of governing equations.
~ Non conservative form
Consider an infinitesimal small fluid element inside nozzle, the fluid same fluid element moving along streamline, such type of condition is non conservative form and the equation obtained in such cases are called as non conservative form of governing equations.
Assumptions for problem:-
- Consider steady isentropic flow through a convergent divergent nozzle.
- Flow is compressible.
- Ignore friction and heat transfer at surface.
- The pressure ratio is maintained to avoid shockwaves.
- The flow is quasi 1D i.e. velocity varies along X-axis only (flow properties are function of x only).
- Length of nozzle = 3.
- The area of nozzle varies according to equation, A=1+2.2(x-1.5)2 0≤x≤3A=1+2.2(x−1.5)2 0≤x≤3
The above figure shows schematic diagram of convergent divergent nozzle.
The flow at inlet of nozzle come from reservoir. For the convergent section flow is subsonic i.e. Mach number < 1, divergent section flow is supersonic i.e. Mach number > 1 and at throat flow is sonic i.e. Mach number = 1.
Where, P0=P0= Pressure at inlet of nozzle.
T0=T0= Temperature at inlet of nozzle.
ρ0=ρ0= Density at inlet of nozzle.
A = Area, V= Velocity, M = Mach number.
P⋆=P⋆= Pressure at throat.
T⋆=T⋆= Temperature at throat.
V⋆=V⋆= a⋆=a⋆= Sound velocity.
Pe=Pe= Pressure at exit of nozzle.
Te=Te= Temperature at exit of nozzle.
Ve=Ve= Velocity at exit of nozzle.
Governing equations:-
As we considering the flow through is only function of X and not Y hence we must apply integral form of Navier-Stokes equation to a control volume consistent with the quasi 1D assumption. Hence by referring below figure we obtain PDE form of governing equations from integral forms.
The shaded region in below figure shows control volume i.e. slice of nozzle flow with infinitesimal thickness of dx.
For the Quasi 1D flow governing equations are as follows-
~Continuity equation
ρ1A1V1=ρ2A2V2ρ1A1V1=ρ2A2V2
~Momentum equation
P1A1+ρ1A1V21+∫A2A1P.da=P2A2+ρ2A2V22P1A1+ρ1A1V21+∫A2A1P.da=P2A2+ρ2A2V22
~Energy equation
h1+V212=h2+V222h1+V212=h2+V222
Where, h=cpTh=cpT
1. Analytical values:-
The mach number can be calculated from below equation.
=> (AA⋆)=1M2[2γ+1(1+γ-12M2)]γ+1γ-1(AA⋆)=1M2[2γ+1(1+γ−12M2)]γ+1γ−1
Hence from above Mach number we can able to calculate pressure, density and temperature values easily.
-> (PP0)=(1+γ-12M2)-γγ-1(PP0)=(1+γ−12M2)−γγ−1
-> (ρρ0)=(1+γ-12M2)-1γ-1(ρρ0)=(1+γ−12M2)−1γ−1
->(TT0)=(1+γ-12M2)-1(TT0)=(1+γ−12M2)−1
By above equations we get exact analytical values for nondimensional pressure ratio, nondimensional density ratio and nondimensional temperature ratio.
The following figures consisting graphs which are plotted by using analytical values v/s nozzle length.
2. Non conservation form:-
~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′]
Where, dash superscript denotes non dimensional variables.
ρ′=ρρ0, A′=AA⋆, T′=TT0, V′=Va0
a0=(γ.R.T0)0.5, t′=tLa0, x′=xL
For computation optimization purpose we have converted Navier-Stokes equation in non dimensional parameters.
3. Conservation form:-
~Continuity equation
∂(ρ′A′)∂t′+∂(ρ′A′V′)∂x′=0
~Momentum equation
∂(ρ′A′V′)∂t′+∂[ρ′A′V′2+(1γ)P′A′]∂x′=(1γ)P′∂A′∂x′
~Energy equation
∂[ρ′(e′γ-1+γ2V′2)A′]dt′+∂[ρ′(e′γ-1+γ2V′2)A′V′+P′A′V′]dx′=0
Let us convert primitive variables in to solution vector 'U', flux vector 'F' and source term 'J'.
Hence we get governing equation as-
∂U1∂t′=-∂F1∂x′
∂U2∂t′=-∂F2∂x′+J2
∂U3∂t′=-∂F3∂x′
Setup of nozzle flow problem:-
Mesh preperation:-
The 1D nozzle domain i.e. X-axis along nozzle can be divided in number of discrete grid points. The last grid which is at nozzle exit can be denote by 'N'. Hence the domain is divided into 'N-1' elements. Here i is arbitrary grid point in domain the grid points i+1 and i+1 are adjacent grid points. The below figure gives clear understanding of mesh preparation. Where Δx denotes mesh element size.
Discretization:-
It can be done by setting a finite difference expressions using MacCormak's explicit technique.
1. Non conservative form:-
Predictor step:[Forward differencing]
~Continuity equation
~Momentum equation
~Energy equation
Now to obtain predicted values of ′ρ′, ′V′ & ′T′
Corrector values:[Backward differencing]
~Continuity equation
~Momentum equation
~Energy equation
The average time derivative values can be calculated by taking average of corrector time derivative and predictor time derivative values as follows,
The final value of ′ρ′, ′V′ & ′T′ can be obtained from following equations.
Initial condition:-
ρ=1-0.3146x
T=1-0.2314x
V=(0.1+1.09x).T12
Boundary conditions:-
->Subsonic inflow condition:-
V1=2V2-V3
->Supersonic outflow condition:-
VN=2VN-1-VN-2
ρN=2ρN-1-ρN-2
TN=2TN-1-TN-2
2. Conservation form:-
Predictor step:
dU1dt=F1(i+1)-F1(i)dx
dU2dt=F2(i+1)-F2(i)dx+J2(i+1)-J2(i)dx
dU3dt=F3(i+1)-F3(i)dx
For corrector method updating values of U1, U2 & U3.
Corrector step:
dU1dt=F1(i)-F1(i-1)dx
dU2dt=F2(i)-F2(i-1)dx+J2(i)-J2(i-1)dx
dU3dt=F3(i)-F3(i-1)dx
The average time derivative values can be calculated by taking average of corrector time derivative and predictor time derivative values as follows,
Hence the final values of U1, U2 & U3 are obtain from following equations.
Initial condition:-
For 0≤x′≤0.5 {ρ′=1T′=1
For 0.5≤x′≤1.5 {ρ′=1-0.366(x′-0.5)T′=1-0.167(x′-0.5)
For 1.5≤x′≤3.5 {ρ′=0.634-0.3879(x′-1.5)T′=0.833-0.3507(x′-1.5)
V′=0.59ρ′A′
Boundary conditions:-
->Subsonic inflow condition:-
U1(i=1)=A′(i=1)
U2(i=1)=2.U2(i=2)-U2(i=3)
U3(i=1)=U1(i=1)(T′i=1γ-1+(γ2)V′2i=1)
->Supersonic outflow condition:-
U1(N)=2U1(N-1)-U1(N-2)
U2(N)=2U2(N-1)-U2(N-2)
U3(N)=2U3(N-1)-U3(N-2)
Timestep calculation (Δt):
By following relation we can able to calculate timestep value.
Δt=CΔxa0+V
Where, C is courant number and by going through some literature review the courant number should be ≤ 1.
Δx is a grid spacing or mesh element size.
a0 is a speed of sound.
So by using some relations the formula for timestep becomes,
Δt′i=CdxT′0.5i+V′i
As there is no certainty for accuracy of local time stepping so for selecting timestep we can use minimum value of timestep in grids.
Δt=minimum(Δtt1,Δtt2,Δtt3,.....,Δtti,......,ΔttN-1,ΔttN)
Although this consistent time marching method may take more time than local time marching method to give converged solution but consistent time marching method gives us physically meaningful transient variation which frequently are of intrinsic value by themselves.
~Non conservation function code:
%-----Non conservative function-----%
function [dt,rho,vel,temp,m_nc,Pressure_nc,mach_nc,nt,rho_t,Pressure_t,m_t,temp_t,vel_t]=N_C(n,x,dx,area,th,c,gamma)
%Initial conditions
rho = 1-0.3146*x;
temp = 1-0.2314*x;
vel = (0.1+1.09*x).*temp.^(0.5);
%Time step calculation
for i=1:n
dt(i)= c*dx/((temp(i))^(0.5)+vel(i));
endfor
dt=min(dt);
nt = 0:1:1400;
%Time loop
for i=1:length(nt)
vel_old=vel;
rho_old=rho;
temp_old=temp;
%Space (grid) loop
%predictor method
for j=2:n-1
dv_dx = (vel(j+1)-vel(j))/dx;
dlnA_dx = (log(area(j+1))-log(area(j)))/dx;
dtemp_dx = (temp(j+1)-temp(j))/dx;
drho_dx = (rho(j+1)-rho(j))/dx;
%Continuity equation
drho_dt_p(j) = -(rho(j)*dv_dx) - (rho(j)*vel(j)*dlnA_dx) - (vel(j)*drho_dx);
%Momentum equation
dv_dt_p(j) = -(vel(j)*dv_dx) - (1/gamma)*(dtemp_dx+(temp(j)/rho(j))*drho_dx);
%Energy equation
dtemp_dt_p(j) = -(vel(j)*dtemp_dx)-(gamma-1)*temp(j)*(dv_dx+vel(j)*dlnA_dx);
%Solution update
rho(j) = rho(j) + drho_dt_p(j)*dt;
vel(j) = vel(j) + dv_dt_p(j)*dt;
temp(j)= temp(j)+ dtemp_dt_p(j)*dt;
endfor
%corrector method
for j=2:n-1
dv_dx = (vel(j)-vel(j-1))/dx;
dlnA_dx = (log(area(j))-log(area(j-1)))/dx;
dtemp_dx = (temp(j)-temp(j-1))/dx;
drho_dx = (rho(j)-rho(j-1))/dx;
%Continuity equation
drho_dt_c(j) = -(rho(j)*dv_dx) - (rho(j)*vel(j)*dlnA_dx) - (vel(j)*drho_dx);
%Momentum equation
dv_dt_c(j) = -(vel(j)*dv_dx) - (1/gamma)*(dtemp_dx+(temp(j)/rho(j))*drho_dx);
%Energy equation
dtemp_dt_c(j) = -(vel(j)*dtemp_dx)-(gamma-1)*temp(j)*(dv_dx+vel(j)*dlnA_dx);
endfor
%Finding average values
drho_dt_nc = 0.5*(drho_dt_p + drho_dt_c);
dv_dt_nc = 0.5*(dv_dt_p + dv_dt_c);
dtemp_dt_nc = 0.5*(dtemp_dt_p + dtemp_dt_c);
%Updating old values
for k=2:n-1
rho(k) = rho_old(k) + drho_dt_nc(k)*dt;
vel(k) = vel_old(k) + dv_dt_nc(k)*dt;
temp(k) = temp_old(k) + dtemp_dt_nc(k)*dt;
endfor
%Boundary condition
%Inflow boundary condition
vel(1) = 2*vel(2)-vel(3);
%Outflow boundary condition
rho(n) = 2*rho(n-1)-rho(n-2);
vel(n) = 2*vel(n-1)-vel(n-2);
temp(n) = 2*temp(n-1)-temp(n-2);
%mass flow rate
m_nc = rho.*area.*vel;
%Pressure
Pressure_nc = rho.*(temp);
%Mach number
mach_nc = vel./(temp).^0.5;
figure(1);
if i==50
plot(x,m_nc,'linewidth',1.75);
hold on;
xlabel('Non dimensional length of nozzle');
ylabel('Non dimensional mass flow rate');
elseif i==100
plot(x,m_nc,'linewidth',1.75);
elseif i==150
plot(x,m_nc,'linewidth',1.75);
elseif i==300
plot(x,m_nc,'linewidth',1.75);
elseif i==700
plot(x,m_nc,'linewidth',1.75);
elseif i==1400
plot(x,m_nc,'linewidth',1.75);
title('Mass flow rate for non conservation form at diffrent timesteps','fontsize',14);
legend('50 timestep','100 timestep','150 timestep','300 timestep','700 timestep','1400 timestep');
end
%Values at throat
rho_t(i)=rho(th);
temp_t(i)=temp(th);
vel_t(i)=vel(th);
Pressure_t(i)=Pressure_nc(th);
m_t(i)=m_nc(th);
endfor
figure(2)
subplot(5,1,1);
plot(nt,rho_t,'linewidth',1.75,'r');
hold on;
xlabel('No of time step');
ylabel('rho / rho_o');
title('Value at throat for variable timesteps (Non conservation form)','fontsize',14);
subplot(5,1,2);
plot(nt,temp_t,'linewidth',1.75,'b');
hold on;
xlabel('No of time step');
ylabel('T / T_o');
subplot(5,1,3);
plot(nt,vel_t,'linewidth',1.75,'g');
hold on;
xlabel('No of time step');
ylabel('Velocity');
subplot(5,1,4);
plot(nt,Pressure_t,'linewidth',1.75,'c');
hold on;
xlabel('No of time step');
ylabel('P / P_o');
subplot(5,1,5);
plot(nt,m_t,'linewidth',1.75,'m');
hold on;
xlabel('No of time step');
ylabel('Massflow rate');
endfunction
All of the above plots shows result for non conservative form.
~Conservation function code:
%-----Conservative function-----%
function [dt,rho_c,vel_c,temp_c,m_c,p_c,mach_c,nt,rho_tc,p_tc,m_tc,temp_tc,vel_tc]=C_f(x,n,dx,area,th,c,gamma)
%Initial conditions
for i=1:length(x)
if x(i)>=0 && x(i)<0.5
rho_c(i) = 1;
temp_c(i) = 1;
elseif x(i)>=0.5 && x(i)<1.5
rho_c(i) = 1-0.366*(x(i)-0.5);
temp_c(i) = 1-0.167*(x(i)-0.5);
elseif x(i)>=1.5 && x(i)<3.5
rho_c(i) = 0.634-0.3879*(x(i)-1.5);
temp_c(i) = 0.833-0.3507*(x(i)-1.5);
endif
endfor
vel_c = (0.59)./(rho_c.*area);
u1 = rho_c.*area;
u2 = rho_c.*area.*vel_c;
u3 = u1.*((temp_c./(gamma-1))+(gamma/2).*vel_c.^2);
%Time step calculation
for i=1:n
dt(i)= (c*dx)/((temp_c(i)^(0.5)) + vel_c(i));
endfor
dt=min(dt);
nt = 0:1:1400;
for i=1:length(nt)
u1_old = u1;
u2_old = u2;
u3_old = u3;
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)));
%predictor step
for j=2:n-1
%Continuity equation
dF1_dx_p(j) = (F1(j+1)-F1(j))/dx;
du1_dt_p(j) = -dF1_dx_p(j);
J2_p(j) = ((gamma-1)/gamma)*(u3(j)-((gamma/2)*(((u2(j))^2)/u1(j))))*((log(area(j+1))-log(area(j)))/dx);
%Momentum equation
dF2_dx_p(j) = (F2(j+1)-F2(j))/dx;
du2_dt_p(j) = -dF2_dx_p(j) + J2_p(j);
%Energy equation
dF3_dx_p(j) = (F3(j+1)-F3(j))/dx;
du3_dt_p(j) = -dF3_dx_p(j);
%updating old values
u1(j) = u1(j) + du1_dt_p(j)*dt;
u2(j) = u2(j) + du2_dt_p(j)*dt;
u3(j) = u3(j) + du3_dt_p(j)*dt;
endfor
rho_c = u1./area;
vel_c = u2./u1;
temp_c = (gamma-1)*((u3./u1)-(gamma/2)*vel_c.^2);
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
%Continuity equation
dF1_dx_c(j) = (F1(j)-F1(j-1))/dx;
du1_dt_c(j) = -dF1_dx_c(j);
J2_c(j) = ((gamma-1)/gamma)*(u3(j)-(gamma/2)*((u2(j))^2/u1(j)))*((log(area(j))-log(area(j-1)))/dx);
%Momentum equation
dF2_dx_c(j) = (F2(j)-F2(j-1))/dx;
du2_dt_c(j) = -dF2_dx_c(j) + J2_c(j);
%Energy equation
dF3_dx_c(j) = (F3(j)-F3(j-1))/dx;
du3_dt_c(j) = -dF3_dx_c(j);
endfor
%Average values
du1_dt = 0.5*(du1_dt_c + du1_dt_p);
du2_dt = 0.5*(du2_dt_c + du2_dt_p);
du3_dt = 0.5*(du3_dt_c + du3_dt_p);
%Updating old values
for k=2:n-1
u1(k) = u1_old(k) + du1_dt(k)*dt;
u2(k) = u2_old(k) + du2_dt(k)*dt;
u3(k) = u3_old(k) + du3_dt(k)*dt;
endfor
%Boundary conditions for problem
%Inflow boundary condition
u1(1)=rho_c(1)*area(1);
u2(1)=2*u2(2) - u2(3);
u3(1)=u1(1)*[((temp_c(1))/(gamma-1))+((gamma/2)*(u2(1)/u1(1))^2)];
%Outflow 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);
rho_c = u1./area;
vel_c = u2./u1;
temp_c = (gamma-1).*((u3./u1)-((gamma/2).*vel_c.^2));
p_c = rho_c.*temp_c;
mach_c = vel_c./((temp_c).^0.5);
m_c = rho_c.*area.*vel_c;
figure(4);
if i==50
plot(x,m_c,'linewidth',1.75);
hold on;
xlabel('Non dimensional length of nozzle');
ylabel('Non dimensional mass flow rate');
elseif i==100
plot(x,m_c,'linewidth',1.75);
elseif i==150
plot(x,m_c,'linewidth',1.75);
elseif i==300
plot(x,m_c,'linewidth',1.75);
elseif i==700
plot(x,m_c,'linewidth',1.75);
elseif i==1400
plot(x,m_c,'linewidth',1.75);
title('Mass flow rate for conservation form at diffrent timesteps','fontsize',14);
legend('50 timestep','100 timestep','150 timestep','300 timestep','700 timestep','1400 timestep');
end
%Values at throat
rho_tc(i) = rho_c(th);
temp_tc(i) = temp_c(th);
vel_tc(i) = vel_c(th);
p_tc(i) = p_c(th);
m_tc(i) = m_c(th);
endfor
figure(5);
subplot(5,1,1);
plot(nt,rho_tc,'linewidth',1.75,'r');
hold on;
xlabel('No of time step');
ylabel('rho/rho_o');
title('Value at throat for variable timesteps (Conservation form)','fontsize',14)
subplot(5,1,2);
plot(nt,temp_tc,'linewidth',1.75,'b');
hold on;
xlabel('No of time step');
ylabel('T / T_o');
subplot(5,1,3);
plot(nt,vel_tc,'linewidth',1.75,'g');
hold on;
xlabel('No of time step');
ylabel('Velocity');
subplot(5,1,4);
plot(nt,p_tc,'linewidth',1.75,'c');
hold on;
xlabel('No of time step');
ylabel('P / P_o');
subplot(5,1,5);
plot(nt,m_tc,'linewidth',1.75,'m');
hold on;
xlabel('No of time step');
ylabel('Massflow rate');
endfunction
All of the above plots shows result for conservative form.
Main octave code:
clear all
close all
clc
%Number of grids
n=31;
%Mesh creation
x=linspace(0,3,n);
dx=x(2)-x(1);
%Cross section area (function of x)
area = 1+2.2*(x-1.5).^2;
%Throat grid calculation
th = find(area==1);
%Courant number
c=0.5;
gamma = 1.4;
gp = gamma+1;
gm = gamma-1;
%radius of nozzle at varius cross section
r=sqrt(area/pi);
%Analytical mach number calculation
for i=1:n
eqn = @(M) (1/M^2)*((2/gp)+((gm*M^2)/gp))^(gp/gm) - area(i)^2;
if ith
x0 = [1+1e-3 5];
M(i) = fzero(eqn,x0);
end
end
%Analytical mach number
mach_ana = M;
%Analytical pressure
Pressure_ana=[1+((gamma-1)/2).*mach_ana.^2].^(-gamma/(gamma-1));
%Analytical density
rho_ana=[1+((gamma-1)/2).*mach_ana.^2].^(-1/(gamma-1));
%Analytical temperature
temp_ana=[1+((gamma-1)/2).*mach_ana.^2].^(-1);
%-----Non conservative function-----%
tic;
[dt,rho,vel,temp,m_nc,Pressure_nc,mach_nc,nt,rho_t,Pressure_t,m_t,temp_t,vel_t]=N_C(n,x,dx,area,th,c,gamma);
ellapsed_nc=toc;
%error for NC form
error_p_nc = abs(Pressure_ana-Pressure_nc);
error_rho_nc = abs(rho_ana-rho);
error_temp_nc = abs(temp_ana-temp);
error_mach_nc = abs(mach_ana-mach_nc);
figure(3);
subplot(4,1,1);
plot(x,Pressure_nc,'linewidth',1.75);
xlabel('Nozzel length');
ylabel('P / P_o');
hold on;
plot(x,Pressure_ana,'linewidth',1.75,'*');
legend('Numerical','Analytical');
title('Non conservation form','fontsize',14)
subplot(4,1,2);
plot(x,rho,'linewidth',1.75);
xlabel('Nozzel length');
ylabel('rho / rho_o');
hold on;
plot(x,rho_ana,'linewidth',1.75,'*');
legend('Numerical','Analytical');
subplot(4,1,3);
plot(x,temp,'linewidth',1.75);
xlabel('Nozzel length');
ylabel('T / T_o');
hold on;
plot(x,temp_ana,'linewidth',1.75,'*');
legend('Numerical','Analytical');
subplot(4,1,4);
plot(x,mach_nc,'linewidth',1.75);
xlabel('Nozzel length');
ylabel('Mach number');
hold on;
plot(x,mach_ana,'linewidth',1.75,'*');
legend('Numerical','Analytical');
%-----Conservative function-----%
tic;
[dt,rho_c,vel_c,temp_c,m_c,p_c,mach_c,nt,rho_tc,p_tc,m_tc,temp_tc,vel_tc]=C_f(x,n,dx,area,th,c,gamma);
ellapsed_c=toc;
%error for Conservation form
error_p_c = abs(Pressure_ana-p_c);
error_rho_c = abs(rho_ana-rho_c);
error_temp_c = abs(temp_ana-temp_c);
error_mach_c = abs(mach_ana-mach_c);
figure(6);
subplot(4,1,1);
plot(x,p_c,'linewidth',1.75);
xlabel('Nozzel length');
ylabel('P / P_o');
hold on;
plot(x,Pressure_ana,'linewidth',1.75,'*');
legend('Numerical','Analytical');
title('Conservation form','fontsize',14)
subplot(4,1,2);
plot(x,rho_c,'linewidth',1.75);
xlabel('Nozzel length');
ylabel('rho / rho_o');
hold on;
plot(x,rho_ana,'linewidth',1.75,'*');
legend('Numerical','Analytical');
subplot(4,1,3);
plot(x,temp_c,'linewidth',1.75);
xlabel('Nozzel length');
ylabel('T / T_o');
hold on;
plot(x,temp_ana,'linewidth',1.75,'*');
legend('Numerical','Analytical');
subplot(4,1,4);
plot(x,mach_c,'linewidth',1.75);
xlabel('Nozzel length');
ylabel('Mach number');
hold on;
plot(x,mach_ana,'linewidth',1.75,'*');
legend('Numerical','Analytical');
%Plot for comparing errors
figure(7);
subplot(4,1,1);
plot(x,error_p_c,'linewidth',1.75);
hold on;
plot(x,error_p_nc,'linewidth',1.75);
xlabel('Nozzel length');
ylabel('pressure error');
legend('Conservation form','Non conservation form');
title('Comparison of error values','fontsize',14);
subplot(4,1,2);
plot(x,error_rho_c,'linewidth',1.75);
hold on;
plot(x,error_rho_nc,'linewidth',1.75);
xlabel('Nozzel length');
ylabel('rho error');
legend('Conservation form','Non conservation form');
subplot(4,1,3);
plot(x,error_temp_c,'linewidth',1.75);
hold on;
plot(x,error_temp_nc,'linewidth',1.75);
xlabel('Nozzel length');
ylabel('Temp error');
legend('Conservation form','Non conservation form');
subplot(4,1,4);
plot(x,error_mach_c,'linewidth',1.75);
hold on;
plot(x,error_mach_nc,'linewidth',1.75);
xlabel('Nozzel length');
ylabel('Mach error');
legend('Conservation form','Non conservation form');
%Plot for comparing non dimensional flow rate
figure(8);
plot(x,m_c,'linewidth',1.75);
hold on;
plot(x,m_nc,'linewidth',1.75);
xlabel('Non dimensional nozzle length');
ylabel('Masss flow rate');
legend('Conservation form','Non-conservation form');
title('Comaprison of mass flow rate','fontsize',14);
%Plot for comparing ellapsed time
figure(9);
bar(1,ellapsed_nc,'r');
hold on;
bar(3, ellapsed_c,'b');
ylabel('Ellapsed time');
title('Comparison for ellapsed time','fontsize',14)
legend('Non-conservative form','Conservative form');
R = max(r);
y =linspace(-R, R, n);
Y = meshgrid(y);
Y = -1*Y';
Mach = zeros(n,n);
for i = 1:n
for j = 1:n
if (Y(j,i) >= -r(i) && Y(j,i) <= r(i))
Mach_nc(j,i) = mach_nc(i);
Mach_c(j,i) = mach_c(i);
else
Mach_nc(j,i) = 0;
Mach_c(j,i) = 0;
end
endfor
endfor
figure(10);
subplot(2,1,1);
colormap(viridis);
contourf(Mach_nc);
title('Non conservative form mach number in nozzle');
subplot(2,1,2);
colormap(viridis);
contourf(Mach_c);
title('Conservative form mach number in nozzle');
For above plot the error for property is calculated by taking difference between analytical value and numerical value. And then that error is plotted on Y-axis against non dimensional nozzle length.
Above plot shows that non dimensional mass flow rate for conservative form shows lots of irregularities than non conservation scheme.
In conservation part we convert our primitive terms into flux terms and then solve it by numerical discretization and again convert those flux values into primitive variables. But in non conservation form we directly solve primitive variables by numerical discretization, hence due to conversion of flux term, source term from primitive variable and after final answer again converting to primitive variables conservation part requires more computational memory and time as compared to the non conservative part for calculations. The above time comparison shows the same i.e. Non conservation form required less time than conservation form.
There are many parameters affect the exact solution and creates error in numerical schemes. Such parameters are listed below-
- Inflow boundary condition error:- At x=0 we have assumed that density, pressure and temperature are equal to reservoir properties, but its strictly possible if Mach number = 0. In above problem as the mach number is function of area ratio so mach number has finite value.
- Courant number:- For linear hyperbolic equation CFL should be ≤0. But as seen in given problem as C increases the error is also increasing. Hence, we have taken C=0.5 for this problem.
- Truncation error/ Grid independence:- Doubling the number of grid points from 31 to 61 changes solution marginally. By comparing the below table, it can be seen that non conservation part provides more close solution to analytic solution than conservation part.
| Flow parameters at nozzle throat | ||||
| ρ⋆ρ0 | T⋆T0 | P⋆P0 | M | |
| Non conservation part n=31 | 0.63866 | 0.83645 | 0.53421 | 0.99939 |
| Conservation part n=31 | 0.69286 | 0.85881 | 0.59504 | 0.91377 |
| Non conservation part n=61 | 0.63766 | 0.83544 | 0.53272 | 0.99977 |
| Conservation part n=61 | 0.63998 | 0.83535 | 0.53461 | 0.99868 |
| Exact analytical solution | 0.63394 | 0.83333 | 0.52828 | 1 |
Conclusion:-
From referring to all above plots and table, we can conclude that for above Quasi 1D nozzle flow problem Non conservation scheme is more accurate, and faster than conservation scheme.
Reference:- A book of 'Computational fluid dynamics' by "John D. Anderson".
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 Smitesh Badnapurkar (41)
Discretization of 3D intake manifold using GEM-3D
Aim:-> To perform discretization of intake manifold with converting it in 1D circuit model from 3D model and observing effect of discretization length on result. Theory: Discretization is the splitting of large parts into smaller sections to improve a model's accuracy. There are two ways in which a fluid…
07 Sep 2021 03:50 PM IST
Conjugate heat transfer analysis of exhaust Port
Aim:-> Conjugate heat transfer (CHT) analysis of exhaust port. Theory:- As we know heat is transferred by 3 modes i.e. Conduction, convection & radiation. Conduction:- is an exchange of energy by direct interaction between molecules of a substance containing temperature differences. Convection:- is…
20 Aug 2021 02:07 PM IST
Turbocharger Modelling using GT Power
Aim:-> Study of turbocharger and its modeling using GT Suit. Theory:- A turbocharger (technically a turbosupercharger), colloquially known as turbo, is a turbine-driven, forced induction device that increases an internal combustion engine's power output by forcing…
20 Aug 2021 08:11 AM IST
Comparison of CI and SI engine and calibration of Single cylinder CI-Engine
Aim:-> To compare SI and CI engine and perform basic calibration of single cylinder CI engine. Comparison of SI V/S CI Engine model: SI Engine CI Engine SI engine is internal combustion engine that operates on the principle of spark ignition. CI engine is internal combustion engine that operates on the principle…
20 Aug 2021 08:09 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.