Week 7 - Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
SIMULATION OF QUASI-ONE-DIMENSIONAL SUPERSONIC NOZZLE FLOW USING MACCORMACK'S TECHNIQUE:
Let us consider the steady, isentropic flow through a convergent-divergent nozzle as show in the following figure:
INLET: Subsonic flow (M<1). Stagnation properties (P0, T0, `rho0`) Area is considered to be infinitely large (A`rarr``infty`) `therefore``V = 0`
THROAT: Sonic flow (M=1)
OUTLET: Supersonic flow (M>1)
THE GOVERNING FLOW EQUATIONS:
Let us consider a control volume,
Net mass flow = outflow - inflow
- (1)
Eqn 1 is the integral form of the continuity equation, where v = Volume, V = Velocity.
The control volume,
p = pressure
`rho`= density
V = velocity
e = internal energy
- (2)
`int int rho Vds` = Outflow - inflow
= `-rhoAV+(rho+drho)(A+dA)(V+dV)`
= `-rhoAV+(rhoA+rhodA+Adrho+dAdrho)(V+dV)`
= `-rhoAV+(rhoAV+rhoVdA+VAdrho+VdAdrho+rhoAdV+rhodAdV+AdrhodV+dAdrhodV))`
`int int rho Vds`= `rhoAdV+rhoVdA+VAdrho`
`int int rho Vds`= `d(rhoAV)` - (3)
From (2) and (3)
`partial/(partialt)(rhoAdx) + d(rhoAV) = 0`
divide by dx,
CONSERVATIVE - DIMENSIONAL FORM:
`(partial(rhoA))/(partialt)+(partial(rhoAV))/(partialx) = 0` `rarr` Continuity equation
`(partial(rhoAV))/(partialt)+(partial(rhoAV^2))/(partialx)=-A(partialp)/(partialx)``rarr` Momentum equation
`(partial(rhoV^2/2A))/(partialt)+(partial(rhoV^3/2A))/(partialx)=-AV(partialp)/(partialx)``rarr` Energy equation
NON CONSERVATIVE - DIMENSIONAL FORM:
`(partial(rhoA))/(partialt)+V(partial(rhoA))/(partialx)=-rhoA(partialV)/(partialx)` `rarr` Continuity equation
`rho(partialV)/(partialt)+rhoV(partialV)/(partialx)=-(partialp)/(partialx)``rarr` Momentum equation
`rho(partiale)/(partialt)+rhoV(partiale)/(partialx)=-p(partialV)/(partialx)-pV(partial(lnA))/(partialx)``rarr` Energy equation
NON-DIMENSIONALIZATION:
MACCORMACK METHOD:
1. CONSERVATIVE - NON DIMENSIONAL FORM:
Continuity equation:
`(partial(rho'A'))/(partialt') + (partial(rho'A'V'))/(partialx')=0`
Momentum equation:
`(partial(rho'A'V'))/(partialt')+(partial[rho'A'V'^2+(1/gamma)P'A'])/(partialx')=1/gammaP'(partialA')/(partialx')`
Energy equation:
`(partial[rho'((e')/(gamma-1)+gamma/2V'^2)A'])/(partialt')+(partial[rho'((e')/(gamma-1)+gamma/2V'^2)A'V'+P'A'V'])/(partialx')=0`
For simplification,
let
`U_1=rho'A'`
`U_2=rho'A'V'`
`U_3=rho'((e')/(gamma-1)+gamma/2V'^2)A'``implies` `U_3=U_1((T')/(gamma-1)+(gamma')/2V'^2)`
`F_1=rho'A'V'`
`F_2=rho'A'V'^2+1/gammaP'A'`
`F_3=rho'((e')/(gamma-1)+gamma/2V'^2)A'V'+P'A'V'`
`J_2=1/gammaP'(partialA')/(partialx')`
The governing equations can be rewritten as:
`(partialU_1)/(partialt')=-(partialF_1)/(partialx')`
`(partialU_2)/(partialt')=-(partialF_2)/(partialx')+J_2`
`(partialU_3)/(partialt')=-(partialF_3)/(partialx')`
writing the values in terms of `U_1,U_2,U_3`
`F_1=U_2`
`F_2=U_2^2/U_1+(gamma-1)/gamma(U_3-gamma/2U_2^2/U_1)`
`F_3=gamma(U_2U_3)/U_1-(gamma(gamma-1))/2U_2^3/U_1^2`
`J_2=(gamma-1)/gamma(U_3-gamma/2U_2^2/U_1)(partial(lnA))/(partialx')`
for simplication the `J_2` term is split into two as follows:
`J_2=J_c(partial(lnA))/(partialx')` where `J_c=(gamma-1)/gamma(U_3-gamma/2U_2^2/U_1)`
Initial conditions:
`V'=U_2/(P'A')=0.59/(rho'A')`
for `1<=x'<=6`
`rho'=1.0`
`T'=1.0`
for `7<=x'<=16`
`rho'=1.0-0.366(x'-0.5)`
`T'=1.0-0.167(x'-0.5)`
for `17<=x'<=31`
`rho'=0.634-0.3879(x'-1.5)`
`T'=0.833-0.3507(x'-1.5)`
Predictor step:
`((partialf)/(partialt))_p` (forward differencing for space terms`
Solution update:
`f^* = f_i+(partialf)/(partialt_p)Deltat`
Corrector step:
`(partialf^*)/(partialt_c)` (backward differencing for space term`
Average:
`((partialf)/(partialt))_(avg) =0.5((partialf)/(partialt_p)+(partialf)/(partialt_c))`
on solving this,
`f^(t+Deltat) = f^t+((partialf)/(partialt))_(avg)Deltat`
The values of `U_1,U_2 &U_3` are calculated using the above method for each time step and the valus of `rho,V&T` are calculated as follows:
`rho'=U_1/(A')`
`V'=U_2/U_1`
`T'=e'=(gamma-1)(U_3/U_1-gamma/2V'^2)`
Boundary conditions:
Inflow:
`U_1(i)=(rho'A')(i)=A'(i)`
`U_2(i)=2U_2(i+1)-U_2(i+2)`
`U_3=U_1((T')/(gamma-1)+(gamma')/2V'^2)`
Outflow:
`U_1(i)=2U_1(i-1)-U_1(i-2)`
`U_2(i)=2U_2(i-1)-U_2(i-2)`
`U_3(i)=2U_3(i-1)-U_3(i-2)`
Separate function is created to calculate `rho`, V, T and m for the conservative - non dimensional form.
function [rho_val v_val t_val m_val u1_val u2_val u3_val rho_t v_t t_t m_t]=conservative_form(n,rho1,v1,t1,u1,u2,u3,a,gamma,dx,dt1,m1)
% Creating a copy
u1_old = u1;
u2_old = u2;
u3_old = u3;
% predictor step
for j=1:n
f1(j)=u2(j);
f2(j)=((u2(j)^2)/u1(j))+(((gamma-1)/gamma)*(u3(j)-((gamma/2)*((u2(j)^2)/u1(j)))));
f3(j)=gamma*((u2(j)*u3(j))/u1(j))-(((gamma*(gamma-1))/2)*((u2(j)^3)/(u1(j)^2)));
jc(j)=((gamma-1)/gamma)*(u3(j)-((gamma/2)*((u2(j)^2)/u1(j))));
endfor
for j = 2:n-1
%forward differencing
du1_dt_p(j)=-(f1(j+1)-f1(j))/dx;
du2_dt_p(j)=(-(f2(j+1)-f2(j))/dx)+(jc(j)*(log(a(j+1))-log(a(j)))/dx);
du3_dt_p(j)=-(f3(j+1)-f3(j))/dx;
%solution update
u1(j)=u1(j)+(du1_dt_p(j))*dt1;
u2(j)=u2(j)+(du2_dt_p(j))*dt1;
u3(j)=u3(j)+(du3_dt_p(j))*dt1;
## rho1(j)=u1_p(j)/a(j);
## v1(j)=u2_p(j)/u1_p(j);
## t1(j)=(gamma-1)*((u3_p(j)/u1_p(j))-((gamma/2)*v1(j)^2));
endfor
%corrector step
for j=2:n-1
f1(j)=u2(j);
f2(j)=((u2(j)^2)/u1(j))+(((gamma-1)/gamma)*(u3(j)-((gamma/2)*((u2(j)^2)/u1(j)))));
f3(j)=gamma*((u2(j)*u3(j))/u1(j))-(((gamma*(gamma-1))/2)*((u2(j)^3)/(u1(j)^2)));
jc(j)=((gamma-1)/gamma)*(u3(j)-((gamma/2)*((u2(j)^2)/u1(j))));
endfor
for j = 2:n-1
%backward differencing
du1_dt_c(j)=-(f1(j)-f1(j-1))/dx;
du2_dt_c(j)=(-(f2(j)-f2(j-1))/dx)+(jc(j)*(log(a(j))-log(a(j-1)))/dx);
du3_dt_c(j)=-(f3(j)-f3(j-1))/dx;
endfor
%computing the average time derivative
for m=2:n-1
du1_dt(m)=0.5*(du1_dt_p(m)+du1_dt_c(m));
du2_dt(m)=0.5*(du2_dt_p(m)+du2_dt_c(m));
du3_dt(m)=0.5*(du3_dt_p(m)+du3_dt_c(m));
%final solution update
u1(m)=u1_old(m)+(du1_dt(m)*dt1);
u2(m)=u2_old(m)+(du2_dt(m)*dt1);
u3(m)=u3_old(m)+(du3_dt(m)*dt1);
% Primitive variables update
rho1(m)=u1(m)/a(m);
v1(m)=u2(m)/u1(m);
t1(m)=(gamma-1)*((u3(m)/u1(m))-((gamma/2)*v1(m)^2));
% mass flow rate
m1(m)=rho1(m)*v1(m)*a(m);
endfor
%boundary conditions
%inflow
u1(1)=rho1(1)*a(1);
u2(1)=(2*u2(2))-u2(3);
v(1) = u2(1)/u1(1);
u3(1)=u1(1)*((t1(1)/(gamma-1))+((v1(1)^2)*(gamma/2)));
## rho1(1)=u1new(1)/a(1);
## v1(1)=u2new(1)/u1new(1);
## t1(1)=(gamma-1)*((u3new(1)/u1new(1))-((gamma/2)*v1(1)^2));
%outflow
u1(31)=(2*u1(30))-u1(29);
u2(31)=(2*u2(30))-u2(29);
u3(31)=(2*u3(30))-u3(29);
## rho1(31)=u1new(31)/a(31);
## v1(31)=u2new(31)/u1new(31);
## t1(31)=(gamma-1)*((u3new(31)/u1new(31))-((gamma/2)*v1(31)^2));
u1_val = u1;
u2_val = u2;
u3_val = u3;
rho_val=rho1;
v_val=v1;
t_val=t1;
m_val=m1;
% variables at throat section
rho_t=rho1(16);
v_t=v1(16);
t_t=t1(16);
% mass flow at throat section
m_t=rho1(16)*v1(16)*a(16);
endfunction2. NON CONSERVATIVE - NON DIMENSIONAL FORM:
Continuity equation:
`(partialrho')/(partialt')=-rho'(partialV')/(partialx')-rho'V'(partial(lnA'))/(partialx')-V'(partialrho')/(partialx')`
Momentum equation:
`(partialV')/(partialt')=-V'(partialV')/(partialx')-1/gamma((partialT')/(partialx')+(T')/(rho')(partialrho')/(partialx'))`
Energy equation:
`(partialT')/(partialt')=-V'(partialT')/(partialx')-(gamma-1)T'[(partialV')/(partialx')+V'(partial(lnA'))/(partialx')]`
Initial conditions:
`rho=1-0.3146x`
`T=1-0.2314x`
`V=(0.1+1.09x)T^(1/2)`
`A=1+2.2(x-1.5)^2` `0<=x<=3`
In general,
Predictor step:
`((partialf)/(partialt))_p` (forward differencing for space terms`
Solution update:
`f^* = f_i+(partialf)/(partialt_p)Deltat`
Corrector step:
`(partialf^*)/(partialt_c)` (backward differencing for space term`
Average:
`((partialf)/(partialt))_(avg) =0.5((partialf)/(partialt_p)+(partialf)/(partialt_c))`
on solving this,
`f^(t+Deltat) = f^t+((partialf)/(partialt))_(avg)Deltat`
Therefor at each grid points the flow field variables are computed as follows.
`rho_i^(t+Deltat) = rho_i^t+((partialrho)/(partialt))_(avg)Deltat`
`V_i^(t+Deltat) = V_i^t+((partialV)/(partialt))_(avg)Deltat`
`T_i^(t+Deltat) = T_i^t+((partialT)/(partialt))_(avg)Deltat`
Time step calculation:
`Deltat=(CFL*dx)/(a+V)`
since a and V varies at each grid point, time step also varies. Therefore the time step is calculated at each grid points and the minimum value is taken as the time step.
`Deltat`=minimum(`Deltat_1^t+Deltat_2^t+....+Deltat_i^t+....+Deltat_n^t``
Boundary conditions:
Inflow:
`V_1=2V_2-V_3`
`rho_1=1`
`T_1=1`
Outflow:
`V_n=2V_(n-1)-V_(n-2)`
`rho_n=2rho_(n-1)-rho_(n-2)`
`T_n=2T_(n-1)-T_(n-2)`
Separate function is created to calculate `rho`, V, T for the non conservative - non dimensional form.
function [rho_val v_val t_val m_val rho_t v_t t_t m_t residual_val1 residual_val2]=non_conservative_form(n,rho2,v2,t2,a,gamma,dx,dt2,m2)
% creating a copy
rho2_old=rho2;
v2_old=v2;
t2_old=t2;
% Predictor step
for j=2:n-1
% forward differencing scheme
drho2dx(j)=(rho2(j+1)-rho2(j))/dx;
dv2dx(j)=(v2(j+1)-v2(j))/dx;
dt2dx(j)=(t2(j+1)-t2(j))/dx;
dlogadx(j)=(log(a(j+1))-log(a(j)))/dx;
% Continuity equation
drho2_dt_p(j)=(-rho2(j)*dv2dx(j))-(rho2(j)*v2(j)*dlogadx(j))-(v2(j)*drho2dx(j));
% Momentum equation
dv2_dt_p(j)=(-v2(j)*dv2dx(j))-((1/gamma)*(dt2dx(j)+((t2(j)/rho2(j))*drho2dx(j))));
% Energy equation
dt2_dt_p(j)=(-v2(j)*dt2dx(j))-((gamma-1)*t2(j)*(dv2dx(j)+(v2(j)*dlogadx(j))));
% Solution update
rho2(j)=rho2_old(j)+drho2_dt_p(j)*dt2;
v2(j)=v2_old(j)+dv2_dt_p(j)*dt2;
t2(j)=t2_old(j)+dt2_dt_p(j)*dt2;
endfor
% Corrector step
for j=2:n-1
% backward differencing scheme
drho2dx(j)=(rho2(j)-rho2(j-1))/dx;
dv2dx(j)=(v2(j)-v2(j-1))/dx;
dt2dx(j)=(t2(j)-t2(j-1))/dx;
dlogadx(j)=(log(a(j))-log(a(j-1)))/dx;
% Continuity equation
drho2_dt_c(j)=(-rho2(j)*dv2dx(j))-(rho2(j)*v2(j)*dlogadx(j))-(v2(j)*drho2dx(j));
% Momentum equation
dv2_dt_c(j)=(-v2(j)*dv2dx(j))-((1/gamma)*(dt2dx(j)+((t2(j)/rho2(j))*drho2dx(j))));
% Energy equation
dt2_dt_c(j)=(-v2(j)*dt2dx(j))-((gamma-1)*t2(j)*(dv2dx(j)+(v2(j)*dlogadx(j))));
endfor
% Computing average time derivative
for m=2:n-1
drho2_dt(m)=0.5*(drho2_dt_p(m)+drho2_dt_c(m));
dv2_dt(m)=0.5*(dv2_dt_p(m)+dv2_dt_c(m));
dt2_dt(m)=0.5*(dt2_dt_p(m)+dt2_dt_c(m));
% final solution update
rho2(m)=rho2_old(m)+(drho2_dt(m)*dt2);
v2(m)=v2_old(m)+(dv2_dt(m)*dt2);
t2(m)=t2_old(m)+(dt2_dt(m)*dt2);
% mass flow rate
m2(m)=rho2(m)*v2(m)*a(m);
endfor
% Applying boundary conditions
% Inlet
v2(1)=2*v2(2)-v2(3);
% Outlet
v2(31)=2*v2(30)-v2(29);
rho2(31)=2*rho2(30)-rho2(29);
t2(31)=2*t2(30)-t2(29);
rho_val=rho2;
v_val=v2;
t_val=t2;
m_val=m2;
% variables at throat section
rho_t=rho2(16);
v_t=v2(16);
t_t=t2(16);
% mass flow at throat section
m_t=rho2(16)*v2(16)*a(16);
% Residual values at throat
residual_val1=drho2_dt(16);
residual_val2=dv2_dt(16);
endfunction
MAIN PROGRAM:
% 1D Quasi Supersonic Nozzle Flow using Maccormack Method
clc
clear all
close all
% domain
L=3;
n=31;
x=linspace(0,L,n);
dx=x(2)-x(1);
gamma=1.4;
% Area profile
a=1+2.2*(x-1.5).^2;
%flow field variables for conservative form (1)
rho1=zeros(1,31);
t1=zeros(1,31);
v1=zeros(1,31);
for i=1:6
rho1(i)=1;
t1(i)=1;
endfor
for i=7:16
rho1(i)=1-(0.366*((x(i)-0.5)));
t1(i)=1-(0.167*((x(i)-0.5)));
endfor
for i=17:31
rho1(i)=0.634-(0.3879*((x(i)-1.5)));
t1(i)=0.833-(0.3507*((x(i)-1.5)));
endfor
for i=1:31
v1(i)=0.59/(rho1(i)*a(i));
endfor
% mass flow rate m1
for i=1:31
m1(i)=rho1(i)*v1(i)*a(i);
end
% Computing Initial conditions
for j = 1:31
u1(j)=rho1(j)*a(j);
u2(j)=u1(j)*v1(j);
u3(j)=u1(j)*((t1(j)/(gamma-1))+((v1(j)^2)*(gamma/2)));
endfor
% Flow-field variables for non conservative form (2)
for i=1:31
rho2(i)=1-(0.3146*x(i)); % Density
t2(i)=1-(0.2314*x(i)); % Temperature
v2(i)=(0.1+1.09*x(i))*(t2(i)^0.5); % Velocity
end
% mass flow rate m2
for i=1:31
m2(i)=rho2(i)*v2(i)*a(i);
end
% Stability criteria
% Calculation of time step based on CFL number
CFL=0.5;
for i=1:n
dt_val1(i)=(CFL*dx)/(((t1(i)).^0.5)+v1(i));
dt_val2(i)=(CFL*dx)/(((t2(i)).^0.5)+v2(i));
endfor
dt1=min(dt_val1);
%dt1 = 0.0267;
dt2=min(dt_val2);
%dt2 = 0.0201;
nt=1400;
time_steps=linspace(1,1400,1400);
% Maccormack method
% Time loop
for k=1:nt
% Conservative form
[rho1_new v1_new t1_new m1_new u1_new u2_new u3_new rho1_t(k) v1_t(k) t1_t(k) m1_t(k)]=conservative_form(n,rho1,v1,t1,u1,u2,u3,a,gamma,dx,dt1,m1);
u1=u1_new;
u2=u2_new;
u3=u3_new;
figure(3)
if k==1
plot(x,m1_new);
end
hold on;
if k==101
plot(x,m1_new);
end
hold on;
if k==201
plot(x,m1_new);
end
hold on;
if k==701
plot(x,m1_new);
end
xlabel('Non dimensional distance through the nozzle - x');
ylabel('Non dimensional mass flow rate');
xticks(0:0.3:3);
legend('0','100','200','700');
% Non - conservative form
[rho2_new v2_new t2_new m2_new rho2_t(k) v2_t(k) t2_t(k) m2_t(k) residual_rho(k) residual_v(k)]=non_conservative_form(n,rho2,v2,t2,a,gamma,dx,dt2,m2);
rho2=rho2_new;
v2=v2_new;
t2=t2_new;
m2=m2_new;
figure(6)
if k==1
plot(x,m2_new);
end
hold on;
if k==51
plot(x,m2_new);
end
hold on;
if k==101
plot(x,m2_new);
end
hold on;
if k==151
plot(x,m2_new);
end
hold on;
if k==201
plot(x,m2_new);
end
hold on;
if k==701
plot(x,m2_new);
end
xlabel('Non dimensional distance through the nozzle - x');
ylabel('Non dimensional mass flow rate');
xticks(0:0.3:3);
legend('0','50','100','150','200','700');
end
% plots for conservative form
figure(1)
subplot(3,1,1);
plot(x,rho1_new);
xlabel('Non dimensional distance through nozzle - x');
ylabel('Non dimensional density');
xticks(0:0.3:3);
title('CONSERVATIVE FORM','color','b');
subplot(3,1,2);
plot(x,v1_new);
xlabel('Non dimensional distance through nozzle - x');
ylabel('Non dimensional velocity');
xticks(0:0.3:3);
subplot(3,1,3);
plot(x,t1_new);
xlabel('Non dimensional distance through nozzle - x');
ylabel('Non dimensional temperature');
xticks(0:0.3:3);
figure(2)
subplot(3,1,1)
plot(time_steps,rho1_t);
xlabel('Time steps');
ylabel('Non dimensional density');
xticks(0:200:1400);
title_text=sprintf('CONSERVATIVE FORM\n Throat condtions');
title(title_text,'color','b');
subplot(3,1,2)
plot(time_steps,v1_t);
xlabel('Time steps');
ylabel('Non dimensional velocity');
xticks(0:200:1400);
subplot(3,1,3)
plot(time_steps,t1_t);
xlabel('Time steps');
ylabel('Non dimensional temperature');
xticks(0:200:1400);
% Non conservative form
figure(4)
subplot(3,1,1)
plot(x,rho2_new,'r');
xlabel('Non dimensional distance through nozzle - x');
ylabel('Non dimensional density');
xticks(0:0.3:3);
title('NON CONSERVATIVE FORM','color','r');
subplot(3,1,2)
plot(x,v2_new,'r');
xlabel('Non dimensional distance through nozzle - x');
ylabel('Non dimensional velocity');
xticks(0:0.3:3);
subplot(3,1,3)
plot(x,t2_new,'r');
xlabel('Non dimensional distance through nozzle - x');
ylabel('Non dimensional temperature');
xticks(0:0.3:3);
figure(5)
subplot(3,1,1)
plot(time_steps,rho2_t,'r');
xlabel('Time steps');
ylabel('Non dimensional density');
xticks(0:200:1400);
title_text=sprintf('NON CONSERVATIVE FORM\n Throat condtions');
title(title_text,'color','b');
subplot(3,1,2)
plot(time_steps,v2_t,'r');
xlabel('Time steps');
ylabel('Non dimensional velocity');
xticks(0:200:1400);
subplot(3,1,3)
plot(time_steps,t2_t,'r');
xlabel('Time steps');
ylabel('Non dimensional temperature');
xticks(0:200:1400);
figure(7)
plot(time_steps,residual_rho,'b');
hold on;
plot(time_steps,residual_v,'r');
xlabel('Time steps');
ylabel('Residual');
xticks(0:200:1400);
legend('Residual density','Residual velocity');
title_text=sprintf('NON CONSERVATIVE\n Residual plot');
title(title_text,'color','b');
% comparison of mass flow in non conservative and conservative forms
figure(8)
plot(x,m1_new,'b');
hold on;
plot(x,m2_new,'r');
xlabel('Non dimensional distance through nozzle - x');
ylabel('Non dimensional mass flow');
xticks(0:0.3:3);
legend('Conservative form','Non conservative form');PLOTS AND DISCUSSIONS:
1. STEADY STATE DISTRIBUTION OF PRIMITIVE VARIABLES THROUGH THE NOZZLE:
- It is observed that in both the forms that density and temperature values decreases along the length of the nozzle from inlet to outlet.
- But the velocity increases from nozzle inlet to outlet.
2. TIME-WISE VARIATION OF THE PRIMITIVE VARIABLES (AT NOZZLE THROAT):
- It is observed that in both the forms there are huge fluctuations in the primitive variables during first few time steps because the flow field variables are driven by a strong potential.
- The steady state is obtained after 700 time steps.
- As the steady state approaches, the changes become smaller.
The variation of the time derivatives with time itself is analysed for non conservative form and the following residual plot is obtained.
- At early times, there is huge fluctuations in the time derivatives because of the unsteady compression and expansion of waves propagating through the nozzle during the transient process.
- At later time steps the change is decreasing and the steady state is obtained after 700 time steps.
3. VARIATION OF MASS FLOW RATE DISTRIBUTION INSIDE THE NOZZLE AT DIFFERENT TIME STEPS:
- After 100 time steps it is observed that there is a huge fluctuation in the massflow rate through the nozzle at different points and starts to settle down around 200 time steps.
- The fluctuations are due to the transient variation of the flow field variables.
- After 700 time steps the mass flow distribution is almost a straight, horizontal line which means that the mass flow has converged to a constant, steady state value throughout the nozzle.
- After 50 time steps it is observed that there is a huge fluctuation in the massflow rate distribution through the nozzle and starts to settle down around 200 time steps.
- The fluctuations are due to the transient variation of the flow field variables.
- After 700 time steps the mass flow distribution is almost a straight, horizontal line which means that the mass flow has converged to a constant, steady state value throughout the nozzle.
4.COMPARISON OF NORMALIZED MASS FLOW RATE DISTRIBUTIONS OF BOTH FORMS:
- In Conservative form, the steady state mass flow rate distribution through the nozzle is almost preserved and remains the same throughout the length of the nozzle. The reason is the mass flow is also a dependent variable in the equations and it is obtained as a primary result.
- However, in Non Conservative form the mass flow rate distribution fluctuates to a greater extent because the mass flow rate is calculated from the primitive variables and hence it is obtained as a secondary result.