Week 7 - Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
Simulation of a 1D Super-sonic nozzle flow using Macormack Method
Introduction
A qausi - 1D numerical analysis is done on a nozzle undergoing subsonic-supersonic isentropic flow. Maccormack’s technique is used for the solving of the discretized governing equation in conservative and non-conservative form.
Convergent- Divergent (C-D) Nozzle – Introduction to physics of the problem
The C-D nozzle has wide variety of application in steam turbines, jet engines and similar applications where axial thrust is needed. The nozzle basically converts the heat energy of a hot pressurized gas into kinetic energy. There are no moving parts involved and this conversion is majorly due to the shape of the nozzle. The thermodynamic and flow properties of the fluid are also responsible for the energy conversion.
A reservoir precedes the nozzle which has very large area compared to it. The hot and high pressure gas at the entry of the nozzle has a negligible velocity. The gas expands isentropically with speeds reaching supersonic levels. The flow is subsonic in the convergent section, increasing to sonic at the throat and gets supersonic at the outlet.
Quasi – one dimensional nozzle flow
The C-D nozzle is at minimum a 2-dimensional flow, i.e the flow changes in x and y directions. For the purpose of this numerical analysis it is assumed that the flow varies only in x dimension. The area of the nozzle can be expressed as a function of x and thus the properties in our assumption will remain constant across the cross section. This type of flow is called as Quasi – one dimensional flow.
Non - conservative form of governing equations
The governing equations can be expressed in the non-dimensional, non-conservative form as follows:-
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 , `T^'=T/T_0 ,ρ^'=ρ/ρ_0 , x^'=x/L ,V^'=V/(√(γRT_0 ) )`
NOTE:- The notations have their usual meanings.
Conservative form of governing equations
Conservative non-dimensional form of equations:-
Continuity equation:-
`(∂U_1)/(∂t^' )= -(∂F_1)/(∂x^' )`
Momentum equation:-
`(∂U_2)/(∂t^' )= -(∂F_2)/(∂x^' ) +J_2`
Energy equation:-
`(∂U_3)/(∂t^' )= -(∂F_3)/(∂x^' )`
Where,
`U_1=ρ^' A'`
`U_2=ρ^' A'V’`
`U_3=ρ^' (T^'/(γ-1)+ γ/2 V'^2 )A'`
`F_1= ρ^' A'V’`
`F_2= ρ^' A^' V'^2+(1/γ ρ^' A^' )`
`F_3=ρ^' (T^'/(γ-1)+γ/2 V'^2 ) V^' A^'+ ρ^' A'V’`
`J_2=1/γ ρ'(∂A')/(∂x^' )`
The equations of conservative form are expressed in solution vector U , the flux vector F and the source term J. The equations are intentionally written in these forms as some instabilities creep in for time-marching solution when written in terms of primitive variables. The flux vector are furthermore expressed in terms of solution vector to negate the numerical instabilities.
`F_1= U_2`
`F_2=(U_2^2)/U_1 +(γ-1)/γ (U_3-γ/2 (U_2^2)/U_1 )`
`F_3=γ (U_2 U_3)/U_1 -γ(γ-1)/2 (U_2^3)/(U_1^2 )`
`J_2=1/γ ρ' T'(∂A')/(∂x^' )`
Boundary conditions for both conservation and non-conservation form
Inflow:-
Velocity is allowed to change to match the mass flow rate at the steady state
The velocity at inflow is extrapolated form the second and third grid point.
`V_1=2V_2-V_3`
The density and temperature are kept constant
`ρ_1=1`
`T_1=1`
Outflow:-
All the variables are allowed to change as follows
`V_N=2V_(N-1)-V_(N-2)`
`ρ_N=2ρ_(N-1)-ρ_(N-2)`
`T_N=2T_(N-1)-T_(N-2)`
Initial conditions for non-conservative form is
`ρ=1-0.3146x`
`T=1-0.2314x`
`V=(0.1+1.09x) T^(1/2)`
Initial conditions for conservative form
Primitive variables for `0≤x^'≤0.5`
`ρ^'=1.0`
`T^'=1.0`
Primitive variables for `0.5≤x^'≤1.5`
`ρ^'=1.0-0.366(x^'-0.5)`
`T^'=1.0-0.366(x^'-0.5)`
Primitive variables for `1.5≤x^'≤3`
`ρ^'=0.634-0.3879(x^'-1.5)`
`T^'=0.833-0.3507(x^'-1.5)`
Calculation of time-step for both conservation and non-conservation form
Courant-Friedrichs-Lowry (CFL) which gives stability criterion for a linear hyperbolic equation can be still used as a guiding equation to determine the timestep for this non-linear partial differential equation problem in a general sense. The Courant number is taken 0.5 to be on the safer side.
`Δ_i^t=(CΔx)/(√(T_i ) +V_i )`
The CFL condition is applied at all grid points for every time step and the minimum time step calculated from the above formula is selected.
Maccormack’s Technique
This technique is used for numerically solving hyperbolic partial differential equation. When a simple upwinding scheme is used it results in numerical dissipation and if higher order accurate schemes are used it may lead to numerical dispersion. This can be avoided by using Maccormack's technique.
A simple algorithm is followed for it as follows:-
TIME LOOP
- Initial conditions
- Predictor step
- Forward difference for space terms:-
`(∂f)/(∂t_p )`
- Solution update:-
Initial conditions + `(∂f)/(∂t_p ) * Δt`
- Corrector step
- Backward difference for space terms:-
`(∂f)/(∂t_c )`
- Average values
`(∂f)/(∂t_(avg) )=1/2 ((∂f)/(∂t_p )+ (∂f)/(∂t_c ))`
- Final update
Initial conditions + `(∂f)/(∂t_(avg) ) * Δt`
END OF TIME LOOP
MAIN CODE
This is the script which calls the conservative and non-conservative form functions and solves them for given value of grid points, gamma and Courant number.
clear all
close all
clc
%inputs
n =61;
gamma = 1.4;
c =0.5;
x = linspace(0,3,n);
dx = x(2)-x(1);
%%non-conservative form
tic
[massflowrate_nc,rho_nc,t_nc,M_nc,p_nc] = non_conservative_form(n,x,dx,gamma,c);
time_nc = toc;
tic
[massflowrate_c,rho_c,t_c,M_c,p_c] = conservative_form(n,x,dx,gamma,c);
time_c = toc;
%comparison of mass flow rates
figure(7);
set(gcf,'Position',[10 10 1000 1000]);
plot(x,massflowrate_nc,x,massflowrate_c);
hold on;
ana = 0.579*ones(1,n);
plot(x,ana,'k');
legend('Non-conservative form','Conservative form','Analytical result','location','north');
grid minor;
xlabel('x/L','fontweight','bold');
ylabel('Normalised mass flow rate','fontweight','bold');
saveas(gcf,'Normalised mass flow rate comparison.png');
FUNCTION FOR NON-CONSERVATION FORM
function[massflowrate,rho_ts,t_ts,M_ts,p] = non_conservative_form(n,x,dx,gamma,c)
rho = 1- 0.3146*x;
t = 1-0.2314*x;
v = (0.1+1.09*x).*t.^0.5;
% area
a = 1+2.2*(x-1.5).^2;
%Time steps
nt = 1400;
%dt = 0.01;
%throat
th = (n+1)/2;
mass_flow_rate = zeros(6,31);
%Outer time loop
for k = 1:nt
%initial conditions
rho_old = rho;
v_old = v;
t_old = t;
%predictor step
for z = 1:n
timestep(z) = c*dx/((t_old(z))^0.5 + v_old(z));
end
dt = min(timestep);
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;
dlogdx = (log(a(j+1))-log(a(j)))/dx;
%continuity equation
drhodt_p(j) = -rho(j)*dvdx - rho(j)*v(j)*dlogdx - 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)*dlogdx);
%solution update
rho(j) = rho_old(j) + drhodt_p(j)*dt;
v(j) = v_old(j)+ dvdt_p(j)*dt;
t(j) = t_old(j) + dtdt_p(j)*dt;
end
%corrector step
for j = 2:n-1
dvdx = (v(j)-v(j-1))/dx;
drhodx = (rho(j)-rho(j-1))/dx;
dtdx = (t(j)-t(j-1))/dx;
dlogdx = (log(a(j))-log(a(j-1)))/dx;
%continuity equation
drhodt_c(j) = -rho(j)*dvdx - rho(j)*v(j)*dlogdx - v(j)*drhodx;
%momentum equation
dvdt_c(j) = -v(j)*dvdx - (1/gamma)*(dtdx + ((t(j)/rho(j))*drhodx));
%energy equation
dtdt_c(j) = -v(j)*(dtdx) - (gamma-1)*t(j)*(dvdx + v(j)*dlogdx);
end
% average time derivative
drhodt = 0.5*(drhodt_c+drhodt_p);
dvdt = 0.5*(dvdt_c+dvdt_p);
dtdt = 0.5*(dtdt_c + dtdt_p);
%final solution update
for i = 2:n-1
v(i) = v_old(i)+ dvdt(i)*dt;
rho(i) = rho_old(i) + drhodt(i)*dt;
t(i) = t_old(i)+ dtdt(i)*dt;
end
%Applying boundary conditions
%inlet
v(1) = 2*v(2) - v(3);
%outlet
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);
% plot(x,v);
% title(sprintf('Time step = %d',k));
% drawnow;
% properties at throat
rho_ts(k) = rho(th);
v_ts(k) = v(th);
t_ts(k) = t(th);
M_ts(k) = v(th)/sqrt(t(th));
p(k) = rho(th)*t(th);
%mass flow rate
if k ==1
mass_flow_rate(1,1:n) =v.*rho.*a;
elseif k == 50
mass_flow_rate(2,1:n) =v.*rho.*a;
elseif k == 100
mass_flow_rate(3,1:n) = v.*rho.*a;
elseif k == 150
mass_flow_rate(4,1:n) = v.*rho.*a;
elseif k == 200
mass_flow_rate(5,1:n) = v.*rho.*a;
elseif k == 700
mass_flow_rate(6,1:n) = v.*rho.*a;
end
end
massflowrate = v.*rho.*a;
%% Timewise variation of variables at nozzle throat
figure(2)
set(gcf,'Position',[10 10 800 800]);
subplot(2,2,1)
plot(rho_ts);
hold on;
plot(1400,0.634,'*');
grid minor;
error1 = ((abs(0.634 - rho_ts(nt)))/0.634)*100;
title({['Normalised Density'],['Error at steady state = ' num2str(error1) '%' ]});
xlabel('Timeteps','fontweight','bold');
ylabel('rho/rho_{0}','fontweight','bold');
subplot(2,2,2)
plot(t_ts);
hold on;
plot(1400,0.833,'*');
error2 = ((abs(0.833 - t_ts(nt)))/0.833)*100;
title({['Normalisaed temperature'],['Error at steady state = ' num2str(error2) '%' ]});
grid minor;
xlabel('Timesteps','fontweight','bold');
ylabel('T/T_0 ','fontweight','bold');
subplot(2,2,3);
plot(M_ts);
hold on;
plot(1400,1,'*');
grid minor;
error3 = ((abs(1 - M_ts(nt)))/1)*100;
title({['Mach number'],['Error at steady state = ' num2str(error3) '%' ]});
xlabel('Timesteps','fontweight','bold');
ylabel('Mach Number','fontweight','bold');
subplot(2,2,4);
plot(p);
hold on;
plot(1400,0.528,'*');
grid minor;
error4 = ((abs(0.528 - p(nt)))/0.528)*100;
title({['Normalised pressure'],['Error at steady state = ' num2str(error4) '%' ]});
xlabel('Timesteps','fontweight','bold');
ylabel('P/P_0 ','fontweight','bold');
saveas(gcf,'Time variation of properties at nozzle throat - non conservative form.png');
%% mass flow rate
figure(3)
set(gcf,'Position',[10 10 800 800]);
plot(x,mass_flow_rate(1,1:n),'b');
hold on;
plot(x,mass_flow_rate(2,1:n),'c');
hold on
plot(x,mass_flow_rate(3,1:n),'*');
hold on;
plot(x,mass_flow_rate(4,1:n),'o');
hold on;
plot(x,mass_flow_rate(5,1:n),'--');
hold on;
plot(x,mass_flow_rate(6,1:n),'k');
hold on;
xlabel('Non-dimensionalised distance through nozzle (x)','fontweight','bold');
ylabel('Non-dimensional Mass flow rate','fontweight','bold');
hold on
plot(3,0.579,'marker','o','markerfacecolor','r');
legend('Deltat 0','Deltat 50','Deltat 100','Deltat 150','Deltat 200','Deltat 700','Analytical result');
grid minor;
saveas(gcf,'Mass flow rate vs time - non conservative form.png');
if n == 31
%% Validation of results
%Validation Data
Rhofinal=[1,0.998,0.997,0.994,0.992,0.987,0.982,0.974,0.963,0.947,0.924,0.892,0.849,0.792,0.721,0.639,0.551,0.465,0.386,0.318,0.262,0.216,0.179,0.150,0.126,0.107,0.092,0.079,0.069,0.061,0.053];
Tfinal=[1,0.999,0.999,0.998,0.997,0.995,0.993,0.989,0.985,0.978,0.969,0.956,0.937,0.911,0.878,0.836,0.789,0.737,0.684,0.633,0.585,0.541,0.502,0.467,0.436,0.408,0.384,0.362,0.342,0.325,0.308];
Vfinal=[0.099,0.112,0.125,0.143,0.162,0.187,0.215,0.251,0.294,0.346,0.409,0.485,0.575,0.678,0.793,0.914,1.037,1.155,1.263,1.361,1.446,1.519,1.582,1.636,1.683,1.723,1.759,1.789,1.817,1.839,1.862];
figure(1);
set(gcf,'Position',[10 10 800 800]);
%SS Plot for rho
subplot(3,1,1)
plot(x,Rhofinal,x,rho,'*');
legend('Validation data','Result from code');
xlabel('x/L','fontweight','bold');
ylabel('rho/rho_0','fontweight','bold');
grid minor;
%SS plot for t
subplot(3,1,2)
plot(x,Tfinal,x,t,'*');
xlabel('x/L','fontweight','bold');
ylabel('T/T_0','fontweight','bold');
legend('Validation data','Result from code');
grid minor;
%SS plot for v
subplot(3,1,3)
plot(x,Vfinal,x,v,'*');
xlabel('x/L','fontweight','bold');
ylabel('V/a_0','fontweight','bold');
legend('Validation data','Result from code','location','southeast');
grid minor;
saveas(gcf,'Validation of results - non conservative form.png');
endFUNCTION FOR CONSERVATION FORM
function[massflowrate,rho_ts,t_ts,M_ts,p] = conservative_form(n,x,dx,gamma,c)
%% calcualte initial profiles
% pressure, temperature
for i = 1:n
if x(i)<= 0.5
rho(i) =1;
t(i) = 1;
elseif 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);
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
% velocity
a = 1+2.2*(x-1.5).^2;
v = 0.59./(rho.*a);
u1 = rho.*a;
u2 = rho.*a.*v;
%massflowrate at O timestep
mass_flow_rate(1,1:n) = u2;
u3 = rho.*a.*((t./(gamma-1)) + (gamma/2)*v.^2);
%throat
th = (n+1)/2;
%timesteps
nt = 1400;
for q = 1 : nt
%% initial conditions
% solution vector
u1_old = u1;
u2_old = u2;
u3_old = u3;
v_old = u2_old./u1_old;
rho_old = u1_old./a;
t_old = (gamma-1)*( (u3_old./u1_old) - (gamma/2)*v_old.^2);
%flux
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);
% % timestep
for z = 1:n
timestep(z) = (c*dx)/((t_old(z))^0.5 + v_old(z));
end
dt = min(timestep);
% dt = 0.0067;
%% predictor step
for j = 2 :n-1
j2(j) = (1/gamma)*rho(j)*t(j)*(a(j+1)-a(j))/dx;
du1dt_p(j) = - (f1(j+1)-f1(j))/dx;
du2dt_p(j) = j2(j)- (f2(j+1)-f2(j))/dx ;
du3dt_p(j) = - (f3(j+1)-f3(j))/dx;
%solution update
u1(j) = u1(j) + du1dt_p(j)*dt;
u2(j) = u2(j) + du2dt_p(j)*dt;
u3(j) = u3(j) + du3dt_p(j)*dt;
end
v_ = u2./u1;
rho= u1./a;
t= (gamma-1)*( (u3./u1) - (gamma/2)*v.^2);
%flux
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 k = 2 :n-1
j2(k) = (1/gamma)*rho(k)*t(k)*((a(k)-a(k-1))/dx);
du1dt_c(k) = - (f1(k)-f1(k-1))/dx;
du2dt_c(k) = j2(k)- (f2(k)-f2(k-1))/dx ;
du3dt_c(k) = - (f3(k)-f3(k-1))/dx;
end
% average values
du1dt_avg = 0.5*(du1dt_p+du1dt_c);
du2dt_avg = 0.5*(du2dt_p+du2dt_c);
du3dt_avg = 0.5*(du3dt_p+du3dt_c);
%final update
for l = 2:n-1
u1(l) = u1_old(l) + du1dt_avg(l)*dt;
u2(l) = u2_old(l) + du2dt_avg(l)*dt;
u3(l) = u3_old(l) + du3dt_avg(l)*dt;
end
% boundary conditions
%inlet
u1(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/2)*v(1)^2);
%outlet
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);
v = u2./u1;
rho = u1./a;
t= (gamma-1)*( (u3./u1) - (gamma/2)*v.^2);
% properties at throat
rho_ts(q) = rho(th);
v_ts(q) = v(th);
t_ts(q) = t(th);
M_ts(q) = v(th)/sqrt(t(th));
p(q) = rho(th)*t(th);
%mass flow rate
if q == 50
mass_flow_rate(2,1:n) =v.*rho.*a;
elseif q == 100
mass_flow_rate(3,1:n) = v.*rho.*a;
elseif q == 150
mass_flow_rate(4,1:n) = v.*rho.*a;
elseif q == 200
mass_flow_rate(5,1:n) = v.*rho.*a;
elseif q == 700
mass_flow_rate(6,1:n) = v.*rho.*a;
end
end
massflowrate = v.*rho.*a;
%mass flow rate vs time
figure(5);
set(gcf,'Position',[10 10 800 800]);
plot(x,mass_flow_rate(1,1:n),'b');
hold on;
plot(x,mass_flow_rate(2,1:n),'c');
hold on
plot(x,mass_flow_rate(3,1:n),'*');
hold on;
plot(x,mass_flow_rate(4,1:n),'o');
hold on;
plot(x,mass_flow_rate(5,1:n),'--');
hold on;
plot(x,mass_flow_rate(6,1:n),'k');
hold on;
xlabel('Non-dimensionalised distance through nozzle (x)','fontweight','bold');
ylabel('Non-dimensional Mass flow rate','fontweight','bold');
hold on
plot(3,0.579,'marker','o','markerfacecolor','r');
legend('Deltat 0','Deltat 50','Deltat 100','Deltat 150','Deltat 200','Deltat 700','Analytical result');
grid minor;
saveas(gcf,'MFR vs time _conser form.png');
%Time variation of properties at throat
figure(6);
set(gcf,'Position',[10 10 800 800]);
subplot(2,2,1);
plot(rho_ts);
hold on;
plot(1400,0.634,'*');
error1 = ((abs(0.634 - rho_ts(nt)))/0.634)*100;
title({['Normalised density'],['Error at steady state = ' num2str(error1) '%' ]});
grid minor;
xlabel('Timeteps','fontweight','bold');
ylabel('rho/rho_{0}','fontweight','bold');
legend('Numerical result', 'analytical result');
subplot(2,2,2)
plot(t_ts);
hold on;
plot(1400,0.833,'*');
error2 = ((abs(0.833 - t_ts(nt)))/0.833)*100;
title({['Normalised temperature'],['Error at steady state = ' num2str(error2) '%' ]});
grid minor;
xlabel('Timesteps','fontweight','bold');
ylabel('T/T_0 ','fontweight','bold');
legend('Numerical result', 'analytical result');
subplot(2,2,3);
plot(M_ts);
hold on;
plot(1400,1,'*');
legend('Numerical result', 'analytical result');
error3 = ((abs(1 - M_ts(nt)))/1)*100;
title({['Mach number'],['Error at steady state = ' num2str(error3) '%' ]});
grid minor;
xlabel('Timesteps','fontweight','bold');
ylabel('Mach Number','fontweight','bold');
subplot(2,2,4);
plot(p);
hold on;
plot(1400,0.528,'*');
error4 = ((abs(0.528 - p(nt)))/0.528)*100;
title({['Normalised pressure'],['Error at steady state = ' num2str(error4) '%' ]});
grid minor;
xlabel('Timesteps','fontweight','bold');
ylabel('P/P_0 ','fontweight','bold');
legend('Numerical result', 'analytical result');
saveas(gcf,'Time variation of properties at throat - conservative form.png');
if n == 31
%% Validation of results
%Validation Data
Rhofinal=[1,0.999,0.997,0.995,0.992,0.988,0.982,0.974,0.962,0.946,0.923,0.891,0.847,0.789,0.718,0.648,0.548,0.462,0.384,0.316,0.260,0.214,0.177,0.148,0.124,0.106,0.090,0.078,0.068,0.060,0.052];
Tfinal=[1,0.999,0.999,0.998,0.997,0.995,0.993,0.989,0.985,0.978,0.969,0.955,0.935,0.909,0.874,0.839,0.783,0.731,0.679,0.628,0.581,0.538,0.500,0.466,0.436,0.409,0.384,0.363,0.344,0.327,0.310];
Vfinal=[0.098,0.110,0.124,0.141,0.161,0.184,0.213,0.249,0.292,0.344,0.408,0.485,0.577,0.682,0.798,0.904,1.046,1.164,1.272,1.368,1.452,1.524,1.586,1.639,1.685,1.725,1.760,1.790,1.817,1.840,1.863];
figure(4);
set(gcf,'Position',[10 10 800 800]);
%SS Plot for rho
subplot(3,1,1)
plot(x,Rhofinal,x,rho,'*');
legend('Validation data','Result from code');
xlabel('x/L','fontweight','bold');
ylabel('rho/rho_0','fontweight','bold');
grid minor;
%SS plot for t
subplot(3,1,2)
plot(x,Tfinal,x,t,'*');
xlabel('x/L','fontweight','bold');
ylabel('T/T_0','fontweight','bold');
legend('Validation data','Result from code');
grid minor;
%SS plot for v
subplot(3,1,3)
plot(x,Vfinal,x,v,'*');
xlabel('x/L','fontweight','bold');
ylabel('V/a_0','fontweight','bold');
legend('Validation data','Result from code','location','southeast');
grid minor;
saveas(gcf,'Validation of results for conservative form.png');
endRESULTS FOR NON-CONSERVATIVE FORM
1) Validation of results
The results of numerical analysis were validated with the data given for 31 grid points and the match is fairly accurate. The validation is done for all the primitive variables across the nozzle length at the 1400th timestep. As expected from a convergent divergent nozzle, the temperature decreases along the length of the nozzle and the velocity increases.
2)Time variation of properties at nozzle throat
The primitive variables at the nozzle throat are plotted against the timesteps. There is initial variation seen in the variables for about 500 timesteps. This variation is attributed to the transient compression and expansion of the waves which travel through the nozzle. The variables thereafter oscillate very little and it can be said that after 1000 timesteps steady state is achieved. The asterisk denotes the analytical value of the variables. From the graphs it can be concluded that the numerical results satisfactorily match the analytical results with the non-conservative approach.
3)Mass flow rate v/s time
The above graph shows the time wise variation of the mass flow rate in the nozzle. Initially at timestep = 0, the mass flow rate varies sinusoidally. This variation is a result of the product of `rho' * a' *v'`. The mass flow rate inside the nozzle tries to stabilise at around 200th timestep from our arbitrary initial condition. At around 700th timestep the mass flow rate becomes constant throughout the nozzle. Thus obeying the continuity equation. Furthermore, the value at which it becomes stable satisfactorily matches the analytical value denoted by the red dot.
RESULTS FOR CONSERVATIVE FORM
1) Validation of results
The results are validated for the conservative approach. Similar to the non-conservative approach the numerical results are in good agreement with the validation data.
2)Time variation of properties at nozzle throat
The error in primitive variables at the throat location when they reach the steady state is more in comparison to the results obtained from the non-conservative approach. The initial oscillation is also captured in the conservative approach but the oscillation is very less in magnitude when compared to the non-conservative approach. This may be due to the different set of initial conditions given to solve the problem.
3)Mass flow rate v/s time
Mass flow rate in conservative approach has a better initial condition as seen from the graph, its is constant throughout the nozzle. This is possible due to the unique way in which the governing equations were formulated. Though the mass flow rate is initially very close to the steady state result, oscillation in its value across the nozzle is still observed. This is due to the wave propogation effect ,as dicussed in the non-conservative approach result.
COMPARISON OF CONSERVATIVE AND NON-CONSERVATIVE FORM
1) Normalised mass flow rate
The above graph compares the normalised mass flow rate for the conservative and non-conservative approach at timestep = 1400. The variation in the values of mass flow rates across the nozzle is observed due to the highly scaled nature of the graph. Practically the values of both the approaches match the analytical result quite satisfactorily.
The key point that can be concluded from this graph is the high oscillation of mass flow rate in the non-conservative approach. This variation is due to the fact that the governing equations are solved for `rho' ,V'` directly and mass flow rate is the 'by-product' we get from the analysis (pun intended). Whereas for the conservative approach the equations are solved primarily for the mass flow rate and thus it gets conserved easliy.
GRID INDEPENDENT TEST
The grid independent test is done to make sure the results from a numerical analysis are not dependent upon the number of grid points used to solve the equations. In the current analysis the code was run for 31 and 61 grid points and its effect on primary variables was observed.
The above table shows values of primitive variables and their percentage error for number of grid points and the approach used. It is observed that non-conservative approach gives higher accuracy compared to the conservative approach for both values of grid points. For all practical purposes we can conclude that 31 grid points is fairly good assumption as the variation in properties is less than 5% for both the approaches.
CONCLUSION
The following conclusion can be drawn from the above analysis:-
- The Maccormack technique can be used successfully to solve hyperbolic partial differential equation of quasi 1-D subsonic supersonic flow, without the need of any artificial viscosity.
- In terms of mass conservation , the conservative approach is superior but in terms of accuracy of primitive variables ,the non-conservative approach is better.
- The solving time is better for non- conservative approach as observed in the below given table. The reason for it can be the extra step of conversion of primitive variables from flux variables required in the conservation approach.
Note:- Time for solving is strictly limited to solving of the equations , all other auxillary (post-processing) commands were commented while recording it.