Simulation of a 1D Super-sonic nozzle flow simulation using Maccormack Method
Objective ::
To solve the isentropic flow of quasi-1D fluid flow through nozzle in conservative and Non-Conservative forms of governing equation using Maccormack method in matlab.
Required to show ::
1. Steady-state distribution of primitive variables inside the nozzle
2. Time-wise variation of the primitive variables
3. Variation of Mass flow rate distribution inside the nozzle at different time steps during the time-marching process
4. Comparison of Normalized mass flow rate distributions of both forms.
Theory ::
Nozzle is used to accelerate or decelerate the fluid flow , to the required range .
1. If the inlet is Subsonic , and outlet is supersonic flow , then the converging section act as a accelerator , diverging section acts as accelerator of flow.
This convergent - Divergent Nozzle , also called as "de laval Nozzle".
2. If the inlet is supersonic and the otlet is supersonic then convergening-diverging duct called as "supersonic diffuser".
Considered the inlet conditions as a reservoir , with finite area (large) Theoretically(A`->oo`) , flow as less (v`approx0`) , Stagnation conditions (`(rho_0 , p_0)`).
Expanded the flow isentropically to the exit conditions , such that Mach number at the throat (=1).
Conditions at Mach = 1 -> Critical (or) Characteristic conditions, flow at M=1 -> Sonic flow.
These variables used in derving the Non-Dimensional Variables.
Flow , locally subsonic in convergent section.
Flow is locally supersonic in divergent section. (`(p_e,rho_e,T_e,M_e)`)
Maximum mass flow rate through the throat arrived at (M=1).
Mach Number >1 at throat is not possible.
Since flow is expanded isentropically , No shock wave phenomena occurs.
If the flow parameters are assumed to change along the one direction -> Quasi 1-D flow , In genral the flow parameters will depend on the both directions.
As the flow parameters changes with area since here area is not constant , need to derive the Differential form of governing equations for quasi 1-D flow , due to approximations.
Below equations are derived in Non-Dimensional form , Using Non dimensional form it is easy to predict the results of the other variables in the form of numbers.
Reduce Number of parameters , used in prepsring the Model of the physical problem.
Non Conservation Form :
Non-Dimensional Variables ::
Conservation Form :
The above form of governing equations are in Conservative Non-Dimensional form.
Defined th variables as a combination of primitive variables , so that our computational power required less , easy to analyse.
U called as solution vectors , since our numerical solution will give numbers directly interms of U1,U2,U3.
F -> Flux as it representes the flux.
Conservative governing equations in terms of dependent variable :
`frac{delU_1}{delt'} = -frac{delF_1}{delx'}`
`frac{delU_2}{delt'} = -frac{delF_2}{delx'} + J_2`
`frac{delU_3}{delt'} = -frac{delF_3}{delx'}`
In conservative form , we must keep in mind that the governing equations , the dependent variables (obtained numbers) are not primitive variables .
Need to decode the dependent variables and update our primitive variables in further steps.
Decoding for primitive variables ,
`rho' = frac{U_1}{A}`
`V' = frac{U_2}{U_1}`
`T'=e'=(gamma-1)*(frac{U_3}{U_1}-frac{gamma}{2}*V'^2)`
Note :
If we represent the Flux terms interms of primitive variables , we observe a blow up in the solution (instabilities) , instead we represent the Flux terms interms of solution vectors.
`F_1 = rho'*A'*V'`
`implies F_1=U_2`
`F_2 = rho'*A'*V'^2 + frac{1}{gamma}*P'*A'`
`implies F_2 = frac{U_2^2}{U_1} + frac{gamma+1}{gamma}*(U_1 - frac{gamma}{2}*frac{U_2^2}{2*U_1})`
`F_3 = gamma*frac{U_2*U_3}{U_1} - frac{gamma*(gamma-1)}{2}*frac{U_2^3}{U_1^2}`
`J_2 = frac{gamma-1}{gamma}*(U_3-frac{gamma}{2}*frac{U_2^2}{U_1})*frac{del(lnA')}{delx'}`
Theoretically implementing the flux terms interms of primitive variables and solution vector no difference , but in Numerical perspective lot of difference interms of stbility of the solution.
we do have an analytical approach , for isentropic quasi -1D flow , in that case using computation power is not nedded . Here our motto is how to approach the solution using MacCormack method.
MacCormack method :
This method is an explicit-finite difference scheme , second order accuracy in both space and time.
This method is modified version of the Lax-Wendroof method , as it is very difficult to used in computation . Undergoes very more calculations , and also in second order accuracy in both space and time.
Consists 3 steps :
1. Predictor step.
2. Corrector step.
3. Averaging step.
In predictor step , we calculate the primitive variables data using (Forward/Backward first order accurate scheme) of initial data.
In corrector step , we update the primitive variables data using (Forward/Backward first order accurate scheme) of predictor step data.
In averaging step , we get the final required primitive variables data. (Averaged of both results).
Keep in mind , we can use forward/Backward schemes for predictor and corrector step , but can't use the same method for both steps.
1.Predictor step :
2. Corrector step :
3. Averaging step :
Similar steps for the conservative method also.
Boundary and initial conditions :
characteristic lines in Nozzle at inflow and outflow of the nozzzle
Implementing the correct boundary and Initial conditions are very imporatant in understanding the solution of the physical problem using computaion.
Non conservative form :
Using method of characteristics , we can say how many Boundary condtions to make float , and hhow many to kept in the dirichilet form.
Velocity is considered as the float boundary condition at inlet.
other in dirichilet form , since only one characteritic line in the direction of flow.
Outflow Boundary conditions :
At outflow, all the characteritic lines , lie outside the domain , so all the boundary conditions are represented in the Neumann(float) Boundary condition form.
Conservative form :
In case of conservative form stability behaviour of the finite-difference formulation of the governing , is slightly more sensitive , therefore it is very better to start with the more
improved initial conditions for the stable solution.
`U_1 = 0.59` -> represents the mass flow rate physically , taken from the result of analytical solution.
By calculating the primitive variables distribution using above equations , we can calculate the dependent source terms using below relations
Time step calculation ::-
Two methods to calcluate the time step ,
For explicit linear hyperbolic equation , CFL <=1 , for the stability.
CFL Condition = (any quantity)*(time-step with which information would travel)/(element step size)
`Deltat` is the time-step value , time taken for the information to propogate fomr one node to other node.
`Deltax` is the value of step size. x(n) - x(n-1). (Diatance b/w the two nodes).
(a+`V'`) velocity with which the characteristic lines propogate.
Here assumed , CFL=0.5.
`1^(st) method` :
calculate the time step at all grid points using the present data for the next iteration.
`2^nd method` :
Take the minimum of all the time step calculated in the above step.
Similar fashion for both Non-conservative and conservative form for the timw-step calculation.
Main code ::
close all
clear all
clc
% Discretizing the domain
x = linspace(0,3,31);
dx = x(2) - x(1);
n = length(x);
area = 1 + 2.2*(x - 1.5).^2;
%throat
th = find(area==1);
% courant Number
c = 0.5;
L = 3;
g = 1.4;
gp = g+1;
gm = g-1;
tic
[rho,vel,temp,pre,mach,mfr,rho_t,temp_t,vel_t,pre_t,mach_t,mfr_t,m_100,m_150] = Non_conservative_form(n,x,dx,area,th,c,g);
a1 = toc;
s_no = (1:31)';
x=x';
area = area';
rho = rho';
vel = vel';
temp = temp';
pre_nc = pre_nc';
mach_Number = mach_nc';
T =table(s_no,x,area,rho,vel,temp,pre_nc,mach_Number);
% Conservative form
N=31;
C = 0.5;
nt = 1400;
gamma = 1.4;
tic
[rho1,temp1,vel1,pre1,mach1,mfr1,rho1_t,temp1_t,vel1_t,pre1_t,Mach1_t,mfr1_t] = Conservative_form_CODE(x,dx,N,C,nt,gamma);
a2 = toc;
s_no = (1:31)';
x=x';
area = area';
rho = rho1';
vel = vel1';
temp = temp1';
pre_nc = pre1';
mach_Number = Mach1';
mass_flow_rate = mfr1';
T2 =table(s_no,x,area,rho,vel,temp,pre_nc,mach_Number,mass_flow_rate);
%%
clc
figure(8)
%a4 = subplot(8,1,1);
plot(x,mfr,'r','LineWidth',1.5)
hold on
plot(x,mfr1,'b','LineWidth',1.5)
legend({'Non-Conservation Form','Conservation Form'})
%h = legend(a4,{'Non-Conservation Form','Conservation Form'});
%set(h,'Location','northeastoutside')
Non-Conservative Form :
function [rho,vel,temp,pre_nc,mach_nc,mfr_nc,rho_t,temp_t,vel_t,pre_t,mach_t,mfr_t,mfr_100,mfr_150] = Non_conservative_form(n,x,dx,area,th,c,g)
%[rho,vel,temp,pre,mach,mfr,rho_t,temp_t,vel_t,pre_t,mach_t,mfr_t] = Non_conservative_form(n,x,dx,area,th,c,g);
% Initial conditions
rho = (1 - 0.3146*x);
temp = (1 - 0.2314*x);
vel = ((0.1 + 1.09*x).*(temp).^(1/2));
n_t= th;
area = area./area(n_t);
nt = 0:1:1400;
%Time loop
for i = 1:length(nt)
vel_old = vel;
rho_old = rho;
temp_old = temp;
for k = 1:n
dt(k) = c*dx/((temp(k))^(0.5) + vel(k));
if k ==n
dt = min(dt);
end
end
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;
drho_dt_p(j) = -(rho(j)*dv_dx) - (rho(j)*vel(j)*dlnA_dx) - (vel(j)*drho_dx);
dv_dt_p(j) = -(vel(j)*dv_dx) - (1/g)*(dtemp_dx+(temp(j)/rho(j))*drho_dx);
dtemp_dt_p(j) = -(vel(j)*dtemp_dx)-(g-1)*temp(j)*(dv_dx+vel(j)*dlnA_dx);
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;
end
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;
drho_dt_c(j) = -(rho(j)*dv_dx) - (rho(j)*vel(j)*dlnA_dx) - (vel(j)*drho_dx);
dv_dt_c(j) = -(vel(j)*dv_dx) - (1/g)*(dtemp_dx+(temp(j)/rho(j))*drho_dx);
dtemp_dt_c(j) = -(vel(j)*dtemp_dx)-(g-1)*temp(j)*(dv_dx+vel(j)*dlnA_dx);
end
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);
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;
end
vel(1) = 2*vel(2) - vel(3);
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
mfr_nc = rho.*area.*vel;
% Pressure
pre_nc = rho.*temp;
% mach number
mach_nc = vel./(temp).^(0.5);
mach_t = vel(th)./(temp(th)).^(0.5);
if i == length(nt)
figure(3)
plot(x,mach_nc)
hold on
end
% Values at throat
rho_t(i) = rho(n_t);
temp_t(i) = temp(n_t);
vel_t(i) = vel(n_t);
pre_t(i) = pre_nc(n_t);
mfr_t(i) = mfr_nc(n_t);
figh1 = figure(1);
figh1.WindowState = 'Maximized';
if i == 50
plot(x,mfr_nc,'linewidth',1.5);
hold on
elseif i==100
mfr_100 = mfr_nc;
plot(x,mfr_nc,'linewidth',1.5);
elseif i==150
mfr_150 = mfr_nc;
plot(x,mfr_nc,'linewidth',1.5);
elseif i==300
plot(x,mfr_nc,'linewidth',1.5);
elseif i==700
plot(x,mfr_nc,'linewidth',1.5,'LineStyle','--');
xlabel('Non dimensional distance through nozzle')
ylabel('Nondimensional mass flow rate')
title('Non dimensional Mass flow rate at varioud time steps')
legend('50','100','150','300','700','1400')
end
end %for
gtext('steady-flow result')
figh2 = figure(2);
figh2.WindowState = 'Maximized';
subplot(5,1,1)
plot(nt,rho_t,'LineWidth',1.5,'color','r')
hold on
ylabel('rho/rho_{0}')
title('Results at throat for various timesteps of Non conservative form')
subplot(5,1,2)
plot(nt,temp_t,'LineWidth',1.5,'color','b')
ylabel('T/T_{o}')
subplot(5,1,3)
plot(nt,vel_t,'linewidth',1.5,'color','g')
ylabel('vel/a_{o}')
subplot(5,1,4)
plot(nt,pre_t,'linewidth',1.5,'color','c')
ylabel('pressure')
subplot(5,1,5)
plot(nt,mfr_t,'linewidth',1.5,'color','m')
ylabel('mass flow rate at throat')
xlabel('Time step value')
figure(3)
% Density
subplot(3,2,1);
plot(x,rho,'r','Linewidth',1.5)
xlabel('Nondimensional Distance through nozzle(x)','fontweight','bold')
ylabel('Nondimensional Density','fontweight','bold')
title('Density Profile in the nozzle')
grid minor;
% Velocity
subplot(3,2,2);
plot(x,vel,'b','Linewidth',1.5)
xlabel('Nondimensional Distance through nozzle(x)','fontweight','bold')
ylabel('Nondimensional Velocity','fontweight','bold')
title('Velocity Profile in the nozzle')
grid minor;
% Temperature
subplot(3,2,3);
plot(x,temp,'g','Linewidth',1.5)
xlabel('Nondimensional Distance through nozzle(x)','fontweight','bold')
ylabel('Nondimensional Temperature','fontweight','bold')
title('Temperature Profile in the nozzle')
grid minor;
% Pressure
subplot(3,2,4);
plot(x,pre_nc,'m','Linewidth',1.5)
xlabel('Nondimensional Distance through nozzle(x)','fontweight','bold')
ylabel('Nondimensional Pressure','fontweight','bold')
title('Pressure Profile in the nozzle')
grid minor;
% Mach Number
subplot(3,2,5);
plot(x,mach_nc,'c','Linewidth',1.5)
xlabel('Nondimensional Distance through nozzle(x)','fontweight','bold')
ylabel('Nondimensional Mach Number','fontweight','bold')
title('Mach Number Profile in the nozzle')
grid minor;
% Mass flow rate
subplot(3,2,6);
plot(x,mfr_nc,'y','Linewidth',1.5)
xlabel('Nondimensional Distance through nozzle(x)','fontweight','bold')
ylabel('Nondimensional Mass flow rate','fontweight','bold')
title('Mass flow rate Profile in the nozzle')
grid minor;
%{
a4 = subplot()
h = legend(a4,{'Non-Conservation Form','Conservation Form','Analytical
solution'})
set(h,'Location','northeastoutside')
%}
Conservative Form :
function [rho_C,T_C,V_C,p_C,Mach_number_C,mfr_C,rho_throat_C,T_throat_C,V_throat_C,p_throat_C,Mach_number_throat_C,mfr_throat_C] = Conservative_form_CODE(x,dx,N,C,nt,gamma)
% Calculating initial conditions for rho and T
for i = 1:N
if x(i)>=0 && x(i)<=0.5
rho_C(i)=1;
T_C(i)=1;
else if x(i)>=0.5 && x(i)<=1.5
rho_C(i)=1-0.366*(x(i)-0.5);
T_C(i)=1-0.167*(x(i)-0.5);
else if x(i)>=1.5 && x(i)<=3.5
rho_C(i)=0.634-0.3879*(x(i)-1.5);
T_C(i)=0.833-0.3507*(x(i)-1.5);
end
end
end
end
% Area and Throat profile
A = 1 + (2.2*((x-1.5).^2));
th = find(A==1);
% Defining initial velocity and pressure
V_C = 0.59./(rho_C.*A);
p_C = rho_C.*T_C;
% Defining solution vectors U1, U2 and U3
U1 = rho_C.*A;
U2 = rho_C.*A.*V_C;
U3 = rho_C.*A.*((T_C./(gamma-1))+((gamma/2).*(V_C.^2)));
CFL = 0.5;
% Speed of sound is given by
% a = (T)^0.5;
% time step
for z=1:N
dt_calc(z) = (CFL*dx)/((T_C(z)).^0.5 + V_C(z));
end
dt = min(dt_calc);
% Number of time steps
nt = 1400;
t_array=1:1400;
% outer time loop
for k = 1:nt
% Saving a copy of old values
U1_old = U1;
U2_old = U2;
U3_old = U3;
% Flux terms in the form of U1, U2, and U3
F1 = U2;
F2 = ((U2.^2)./(U1)) + ((gamma-1)/gamma).*(U3-(0.5*gamma.*U2.^2)./U1);
F3 = ((gamma.*U2.*U3)./U1)-(0.5*gamma*(gamma-1).*U2.^3)./(U1.^2);
% Predictor step
% Forward differences
for j=2:N-1
dF1_dx = (F1(j+1)-F1(j))/dx;
dF2_dx = (F2(j+1)-F2(j))/dx;
dF3_dx = (F3(j+1)-F3(j))/dx;
% Continuity equation
dU1_dt_p(j) = -dF1_dx;
% Momentum equation
% Source term in the momentum equation
J2(j) = (1/gamma).*rho_C(j).*T_C(j).*((A(j+1)-A(j))/dx);
dU2_dt_p(j) = -dF2_dx + J2(j);
% Energy equation
dU3_dt_p(j) = -dF3_dx;
% Solution update
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;
end
% Updating the values
rho_C = U1./A;
V_C = U2./U1;
T_C= (gamma-1).*((U3./U1)-((gamma./2).*(V_C.^2)));
p_C = rho_C.*T_C;
% Updating the flux terms after the predictor step
F1 = U2;
F2 = ((U2.^2)./(U1)) + ((gamma-1)/gamma).*(U3-(0.5*gamma.*U2.^2)./U1);
F3 = ((gamma.*U2.*U3)./U1)-(0.5*gamma*(gamma-1).*U2.^3)./(U1.^2);
% Corrector step
% Backward differences
for j=2:N-1
dF1_dx = (F1(j)-F1(j-1))/dx;
dF2_dx = (F2(j)-F2(j-1))/dx;
dF3_dx = (F3(j)-F3(j-1))/dx;
% Continuity equation
dU1_dt_c(j) = -dF1_dx;
% Momentum equation
% Source term in the momentum equation
J2(j) = (1/gamma).*rho_C(j).*T_C(j).*((A(j)-A(j-1))/dx);
dU2_dt_c(j) = -dF2_dx + J2(j);
% Energy equation
dU3_dt_c(j) = -dF3_dx;
end
% Computing the average time derivative
dU1_dt_avg = 0.5.*(dU1_dt_p + dU1_dt_c);
dU2_dt_avg = 0.5.*(dU2_dt_p + dU2_dt_c);
dU3_dt_avg = 0.5.*(dU3_dt_p + dU3_dt_c);
% Final solution update
for i = 2:N-1
U1(i) = U1_old(i) + dU1_dt_avg(i)*dt;
U2(i) = U2_old(i) + dU2_dt_avg(i)*dt;
U3(i) = U3_old(i) + dU3_dt_avg(i)*dt;
end
% Apply the boundary conditions
% At Inlet
U1(1) = rho_C(1)*A(1);
U2(1) = 2*U2(2) - U2(3);
V_C(1) = U2(1)./U1(1);
U3(1) = U1(1)*((T_C(1)/(gamma-1)) + ((gamma/2)*(V_C(1)).^2));
% At 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);
% Updation of flow field variables
rho_C = U1./A;
V_C = U2./U1;
T_C = (gamma-1).*((U3./U1)-((gamma./2).*(V_C.^2)));
p_C = rho_C.*T_C;
% Calculating Mach number and mass flow rate (mfr)
Mach_number_C = V_C./(T_C.^0.5);
mfr_C = rho_C.*V_C.*A;
% Values at throat for each iterations
rho_throat_C(k) = rho_C(th);
V_throat_C(k) = V_C(th);
T_throat_C(k) = T_C(th);
p_throat_C(k) = p_C(th);
Mach_number_throat_C(k) = Mach_number_C(th);
mfr_throat_C(k) = mfr_C(th);
% Instantaneous distributions of the nondimensional mass flow through the nozzle
figure(5);
hold on;
if k==50
plot(x,mfr_C,'r','linewidth',1.25);
else if k==100
plot(x,mfr_C,'m','linewidth',1.25);
else if k==150
plot(x,mfr_C,'b','linewidth',1.25);
else if k==200
plot(x,mfr_C,'g','linewidth',1.25);
else if k==700
plot(x,mfr_C,'k','linewidth',1.25);
end
end
end
end
end
end
title('Conservative form - Instantaneous distributions of the nondimensional mass flow through the nozzle');
xlabel('Nondimensional distance through nozzle (x)','fontweight','bold');
ylabel('Nondimensional mass flow (rhoVA/rho_{o}V_{o}A_{o}','fontweight','bold');
legend('50 Deltat','100 Deltat','150 Deltat','200 Deltat','700 Deltat');
axis([0 3 0 2]);
grid on;
grid minor;
% plotting
% Steady-state distributions of nondimensional variables through the nozzle
figure(6)
% Density
subplot(3,2,1);
plot(x,rho_C,'r','Linewidth',1.5)
xlabel('Nondimensional Distance through nozzle(x)','fontweight','bold')
ylabel('Nondimensional Density','fontweight','bold')
title('Density Profile in the nozzle')
grid minor;
% Velocity
subplot(3,2,2);
plot(x,V_C,'g','Linewidth',1.5)
xlabel('Nondimensional Distance through nozzle(x)','fontweight','bold')
ylabel('Nondimensional Velocity','fontweight','bold')
title('Velocity Profile in the nozzle')
grid minor;
% Temperature
subplot(3,2,3);
plot(x,T_C,'b','Linewidth',1.5)
xlabel('Nondimensional Distance through nozzle(x)','fontweight','bold')
ylabel('Nondimensional Temperature','fontweight','bold')
title('Temperature Profile in the nozzle')
grid minor;
% Pressure
subplot(3,2,4);
plot(x,p_C,'c','Linewidth',1.5)
xlabel('Nondimensional Distance through nozzle(x)','fontweight','bold')
ylabel('Nondimensional Pressure','fontweight','bold')
title('Pressure Profile in the nozzle')
grid minor;
% Mach Number
subplot(3,2,5);
plot(x,Mach_number_C,'m','Linewidth',1.5)
xlabel('Nondimensional Distance through nozzle(x)','fontweight','bold')
ylabel('Nondimensional Mach Number','fontweight','bold')
title('Mach Number Profile in the nozzle')
grid minor;
% Mass flow rate
subplot(3,2,6);
plot(x,mfr_C,'r','Linewidth',1.5)
xlabel('Nondimensional Distance through nozzle(x)','fontweight','bold')
ylabel('Nondimensional Mass flow rate','fontweight','bold')
title('Mass flow rate Profile in the nozzle')
grid minor;
% Timewise variations of the flow-field variables at the nozzle throat(where A=1)
figh2 = figure(7);
figh2.WindowState = 'Maximized';
subplot(6,1,1);
plot(t_array,rho_throat_C,'r','LineWidth',1.5)
ylabel('rho/rho_{0}','fontweight','bold')
title('Results at throat for various timesteps of Non conservative form')
subplot(6,1,2);
plot(t_array,V_throat_C,'g','LineWidth',1.5)
ylabel('vel/a_{o}','fontweight','bold')
grid minor
subplot(6,1,3);
plot(t_array,T_throat_C,'b','LineWidth',1.5)
ylabel('T/T_{o}','fontweight','bold')
grid minor
subplot(6,1,4);
plot(t_array,p_throat_C,'c','LineWidth',1.5)
ylabel('Pressure','fontweight','bold')
grid minor
subplot(6,1,5);
plot(t_array,Mach_number_throat_C,'m','Linewidth',1.5)
ylabel('Mach number','fontweight','bold')
grid minor
subplot(6,1,6);
plot(t_array,mfr_throat_C,'r')
xlabel('Number of time steps','fontweight','bold')
ylabel('Mass flow rate','fontweight','bold')
grid minor
xlabel('Time step','fontweight','bold')
end
Results ::
Mass flow rate for various time steps in Non-Conservation-Form :
Observed the steady state at the timestep -> 700.
Non dimensional variables at throat for timesteps in Non Conservation Form:
Non dimensional parameters throughout nozzle at steady state in Non-Conservation form:
Mass flow rate for various time steps throughout nozzle for various timesteps in Conservation form:
Non dimensional variables at throat for timesteps in Conservation Form:
Non dimensional parameters throughout nozzle at steady state in Conservation form:
Mass flow Rate Comparision for both Conservative and Non-Conservative Form :
Can observe from the above plots , Non conservative form are under fluctuations , but in very minute ranges , In conservative form more accurate result than Non-Conservative form
If considered for more accurate solution -> vonservatibe form is preferable.
But for the conservative form , more computation power required , In reality there would be millions of grid cells.
Non-conservative form is robust than conservative , instabilities occur in the Conservative form.
So we need to choose the method carefully with respect to our requirements.
Non-conservative form results in faster convergence compared to the conservative form.
Grid Independence test ::
As increased to the n = 91 grids , solution became unstable , so, it is not good to make a fine mesh than n=61 grids.