Simulation of a 1D Super-sonic nozzle flow using Macormack Method
Introduction
The project concerns with steady isentropic flow through a convergent divergent nozzle. The entry to the nozzle is linked to a reservoir, which is the source of the fluid, concerned. The flow inside the nozzle expands adiabatically from subsonic (M<1) at the convergent zone to sonic (M=1) at the throat and finally to supersonic (M>1) at the divergent section. The illustration of the problem is given below:
From the above figure, the parameters menioned are identified as:
``P_0`= Stagnation Pressure at the inlet of Nozzle
`T_0`= Stagnation Temperature at inlet of Nozzl`
`rho_0`= Density of flow at nozzle inlet
V = Velocity of flow at Nozzle inlet
A = Area of Reservoir at Nozzle inlet
A* = Cross-Secion Area of Nozzle Throat
p* = Pressure of flow at Nozzle Throat
T* = Temperature of flow at Nozzle Throat
V* = Velocity of flow at Nozzle throat at sonic speed (M=1)=`a_0`
M = Mach No. at Nozzle Throat = 1
`p_e`= Exit Pressure of Nozzle Flow
`T_e`= Temperature at exit of Nozzle Flow
`V_e`= Velocity of flow at exit of Nozzle
`M_e`= Mach no. at Nozzle exit
L = Length of the Nozzle
Assumptions:
a. The flow through nozzle is steady.
b. The flow through nozzle is isentropic
c. The area of the reservoir to which the nozzle is connected is very large and tends to infinity.
d. As a corollary of assumption (c), velocity at nozzle entry is equal to nil.
Objectives
1. Numerical Solution of governing equations in both conservative and non conservative form, to calculate flow parameters.
2. Demonstrating time-wise variation of flow parameters at nozzle throat.
3. Demonstrating spacewise variation of flow parameters throughout the nozzle after the completion of the observed time interval
4. Comparison of Non Dimensional mass flow rate between conservative and non conservative form.
5. Performing grid independence test.
6. Computation time for conservative and non conservative methods.
Non Dimensional Variables
Temperature: `T'=T/T_0`
Density: `rho'=rho/rho_0`
Displacement: `x'=x/L`
Velocty: `V'=V/a_0`
Time: `t'=t/(L/a_0)`
Area: `A'=A/(A*)`
The governing equations are expressed in the form of the non dimensional variables, as explained above. The conservative and non conservative forms of the governing equaions, using the non dimensional variables are listed below.
Governing Equations: Non Conservative Form
Continuity Equation:
`(delrho')/(delt')=rho'(delV')/(delx')-rho'V'(dellnA')/(delx')-V'(delrho')/(delx')`
Momentum Equation:
`(delV')/(delt')=V'(delV')/(delx')-(1/gamma)((delT')/(delx')+(T')/(rho')(delrho')/(delx'))`
Energy Equation:
`(delT')/(delt')=V'(delT')/(delx')-(gamma-1)T'((delV')/(delx')+V'(dellnA')/(delx'))`
Governing Equations: Conservative Form
Continuity Equation:
`(del(rho'A'))/(delt')+(del(rho'A'V'))/(delx')=0`
Momentum Equation:
`(del(rho'A'V'))/(delt')+(del(rhoAV'^2+1/(gamma)p'A'))/(delx')=(1/gamma)p'(delA')/(delx')`
Energy Equation:
`(del[rho'(T/(gamma-1)+gamma/2V'^2)A'])/(delt')+(del[rho'(T/(gamma-1)+gamma/2V'^2)V'A'+p'A'V'])/(delx')=0`
Solution Algorithm: MacCormack Method Description
Step 1: From initial condition, value of function at all the grid points are obatined at t=0
Step2: The continuity, momentum and energy equations are solved by forward difference method to obtain `((delf)/(delt))_p`, where f stands for density, velocity of temperature. The p mentioned in this step indicates the fact, that the partial derivatives obtained are in the predictor step.
Step3: Solution update is done for all 3 of density, temperature and velocity by the equation `f*=f+((delf)/(delt))_p*dt`
Step4: The continuity, momentum and energy equations are again solved, this time by backward difference method, using the updated values f*(density, temperature and velocity obtained from Step 3), to obtain `((delf*)/(delt))_c`. The c mentioned in this step indicates the fact, that the partial derivatives obtaind are in the corrector step.
Step5: The mean partial derivative is found out by `((delf)/(delt))_(avg)=0.5(((delf)/(delt))_p+((delf*)/(delt))_c)`
Step6: Finally the flow parameters are obtained at time step 1 (flow parameters at time step 0 were already given as part of the initial conditions in Step 1) by `f_1=f_0+((delf)/(delt))_(avg)*dt` where f_0 denotes the flow properties at time = 0 and is same as f, that is considered in the above 5 steps.
Solution Strategy for Non Conservative Governing Equations
For non conservative equations, the solution strategy is pretty straight forward. The time derivatives are expressed in terms of the primitive variables i.e. density, temperature and velocity. Hence, it is trivial to obtain `(delf)/(delt)`by MacCormack method as illustrated above.
Solution Strategy for Conservative Governing Equations
In conservative forms, the time derivatives are not given in the form of primitive variables. Hence, it is essential to solve the time derivatives first and then solve for the time derivatives of the primitive variables sequentially. Hence, the following substitutions are done
Solution terms:
`U_1=rho'A'`
`U_2=rho'A'V'`
`rho'(T/(gamma-1)+gamma/2V'^2)A'`
Flux Terms:
`F_1=rho'A'V'=U_2`
`F_2=rho'A'V'^2+(1/gamma)p'A'=U_2^2/U_1+(gamma-1)/(gamma)(U_3-(gamma/2)*U_2^2/U_1)`
`F_3=rho'((T')/(gamma-1)+gamma/2V'^2)V'A"+p'A'V'=gamma*(U_2U_3)/U_1-gamma((gamma-1)/2)U_2^3/U_1^2`
Source Term:
`J_2=(1/gamma)p'A'=(gamma-1)/(gamma)(U_3-(gamma/2)*U_2^2/U_1)(del(lnA')/(delx'))`
Primitive Variables expressed as Solution terms:
`rho'=U_1/A`
`V'=U_2/U_1`
`T'=(gamma-1)(U_3/U_1-(gamma/2)V'^2)`
`p'=rho'T'`
Hence, in the conservative form of the equation, the solution is initially done for the solutions terms Ui by MacCormack method and then the primitve variables are solved from the solution terms.
MATLAB Code:
Main Script
clear all
close all
clc
%Input
L=3;
n=31; %No. of grid points
x=linspace(0,L,n); %Domain Length
dx=x(2)-x(1); %Grid spacing
gamma=1.4;
nt=1400;
throat=(n-1)/2;
a=1+2.2*(x-1.5).^2; %Area Profile which remains constant with time
C=0.5;
tic
[m, M,v, Temp, p, rho, Temp_throat, rho_throat, p_throat, M_throat]=non_conservative_form(x, dx, n, C, a, nt, gamma, throat);
time_elapsed=toc;
fprintf('The Computation time by Non-Conservative form: %0.4f',time_elapsed)
%% Time-wise variation of density, temperature, pressure and Mach number for the Non-Conservative form:
time_variation = figure(2);
subplot(4,1,1);
plot(rho_throat,'b');
title('Time-wise variations in quantity in Non-Conservative form');
legend('Density');
ylabel('rho/rho_{o}');
axis([0 1400 0 1.3]);
grid minor;
subplot(4,1,2);
plot(Temp_throat,'r');
legend('Temperature');
ylabel('T/T_{o}');
axis([0 1400 0.6 1]);
grid minor;
subplot(4,1,3);
plot(p_throat,'g');
legend('Pressure');
ylabel('p/p_{o}');
axis([0 1400 0.4 1.1]);
grid minor;
subplot(4,1,4);
plot(M_throat,'m');
legend('Mach Number');
ylabel('M');
axis([0 1400 0.6 1.4]);
grid minor;
xlabel('Number of iterations');
%% Qualitative Aspects of 1D flow in Non-Conservative Form:
qualitative_non_cons = figure(3);
subplot(4,1,1);
plot(x,M,'-b');
title('Qualitative Aspects of 1D Supersonic Nozzle flow: Non-Conservative form');
ylabel('M');
legend('Mach Number','location','southeast');
grid minor;
subplot(4,1,2);
plot(x,p,'-r');
ylabel('p/p_{o}');
legend('Pressure ratio');
grid minor;
subplot(4,1,3);
plot(x,rho,'-g');
ylabel('rho/rho_{o}');
legend('Density ratio','outside');
grid minor;
subplot(4,1,4);
plot(x,Temp,'-m');
ylabel('T/T_{o}');
legend('Temperature ratio','outside');
grid minor;
%% Solution by Conservative form:
tic
[m_c, M_c, v_c, Temp_c, p_c, rho_c, Temp_throat_c, rho_throat_c, p_throat_c, M_throat_c] = conservative_form(x,dx,n,C,a,nt,gamma,throat);
time_elapsed_cons = toc;
fprintf('The computation time by Conservative form: %0.4f',time_elapsed_cons);
%% Time-wise variation of density, temperature, pressure and Mach number for the Non-Conservative form:
time_variation_cons = figure(5);
subplot(4,1,1);
plot(rho_throat_c,'b');
title('Time-wise variations in quantity in Conservative Form');
legend('Density');
ylabel('rho/rho_{o}');
axis([0 1400 0 1.3]);
grid minor;
subplot(4,1,2);
plot(Temp_throat_c,'r');
legend('Temperature');
ylabel('T/T_{o}');
axis([0 1400 0.6 1]);
grid minor;
subplot(4,1,3);
plot(p_throat_c,'g');
legend('Pressure');
ylabel('p/p_{o}');
axis([0 1400 0.4 1.1]);
grid minor;
subplot(4,1,4);
plot(M_throat_c,'m');
legend('Mach Number');
ylabel('M');
axis([0 1400 0.6 1.4]);
grid minor;
xlabel('Number of iterations');
%% Qualitative Aspects of 1D flow in Conservative Form:
qualitative_cons = figure(6);
subplot(4,1,1);
plot(x,M_c,'-b');
title('Qualitative Aspects of 1D Supersonic Nozzle flow: Conservative form');
ylabel('M');
legend('Mach Number','location','southeast');
grid minor;
subplot(4,1,2);
plot(x,p_c,'-r');
ylabel('p/p_{o}');
legend('Pressure ratio');
grid minor;
subplot(4,1,3);
plot(x,rho_c,'-g');
ylabel('rho/rho_{o}');
legend('Density ratio','outside');
grid minor;
subplot(4,1,4);
plot(x,Temp_c,'-m');
ylabel('T/T_{o}');
legend('Temperature ratio','outside');
grid minor;
%% Comparison of Mass-Flow rates by Conservative and Non-Conservative methods and comparison with Exact Form:
mass_comparison = figure(7);
hold on;
plot(x, m, 'g','LineWidth', 1.75);
plot(x, m_c, 'm','LineWidth', 1.75);
line([0 3], [0.579 0.579],'LineWidth', 1.75);
legend('Non-Conservative Form', 'Conservative Form', 'Exact Solution');
grid minor;
xlabel('x/L');
ylabel('rhoVA/rho_{o}V_{o}A_{0}');
title('Non-Dimensional Mass Flow comparison with Exact Solution');
Function to Solve Non Conservative form of Equation
function[m, M,v, Temp, p, rho, Temp_throat, rho_throat, p_throat, M_throat]=non_conservative_form(x, dx, n, C, a, nt, gamma, throat)
rho = 1-0.3146*x; %initial density condition
Temp=1-0.2314*x; %initial temperature condition
v=(0.1+1.09*x).*Temp.^0.5; %initial velocity condition
for k=1:nt
rho_old=rho;
v_old=v;
Temp_old=Temp ;
dt=min(C.*dx./(sqrt(Temp)+v));
%Predictor Step
for j=2:n-1
dvdx=(v(j+1)-v(j))/dx;
drhodx=(rho(j+1)-rho(j))/dx;
dTempdx=(Temp(j+1)-Temp(j))/dx;
dlogadx=(log(a(j+1))-log(a(j)))/dx;
%Continuity Equation
drhodt_p(j)=-rho(j)*dvdx-rho(j)*v(j)*dlogadx-v(j)*drhodx;
%Momentum Equation
dvdt_p(j)=-v(j)*dvdx-(1/(gamma))*(dTempdx+(Temp(j)/rho(j))*drhodx);
%Energy Equation
dTempdt_p(j)=-v(j)*dTempdx-(gamma-1)*Temp(j)*(dvdx+v(j)*dlogadx);
%Solution Update
rho(j)=rho(j)+drhodt_p(j)*dt;
v(j)=v(j)+dvdt_p(j)*dt;
Temp(j)=Temp(j)+dTempdt_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;
dTempdx=(Temp(j)-Temp(j-1))/dx;
dlogadx=(log(a(j))-log(a(j-1)))/dx;
%Continuity Equation
drhodt_c(j)=-rho(j)*dvdx-rho(j)*v(j)*dlogadx-v(j)*drhodx;
%Momentum Equation
dvdt_c(j)=-v(j)*dvdx-(1/(gamma))*(dTempdx+(Temp(j)/rho(j))*drhodx);
%Energy Equation
dTempdt_c(j)=-v(j)*dTempdx-(gamma-1)*Temp(j)*(dvdx+v(j)*dlogadx);
end
drhodt=0.5*(drhodt_c+drhodt_p);
dvdt=0.5*(dvdt_c+dvdt_p);
dTempdt=0.5*(dTempdt_c+dTempdt_p);
for i=2:n-1;
rho(i)=rho_old(i)+drhodt(i)*dt;
v(i)=v_old(i)+dvdt(i)*dt;
Temp(i)=Temp_old(i)+dTempdt(i)*dt;
end
%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);
Temp(n)=2*Temp(n-1)-Temp(n-2);
p=rho.*Temp;
m=rho.*a.*v;
M=v./sqrt(Temp);
Temp_throat(k)=Temp(throat);
p_throat(k)=p(throat);
rho_throat(k)=rho(throat);
M_throat(k)=M(throat);
mass=figure(1);
hold on;
if k == 50
plot(x,m,'b','linewidth',1.75);
else if k == 100
plot(x,m,'g','linewidth',1.75);
else if k == 150
plot(x,m,'r','linewidth',1.75);
else if k == 200
plot(x,m,'m','linewidth',1.75);
else if k == 700
plot(x,m,'-','linewidth',1.75);
plot(x(end),m(end),'b*');
end
end
end
end
end
end
title('Comparison of non-dimensional mass flow rates at different time-instants');
xlabel('x/L (Non-dimensional distance through the nozzle)');
ylabel('rhoVA/rho_{o}V_{o}A_{o}');
legend('50Deltat','100Deltat','150Deltat','200Deltat','700Deltat');
axis([0 3 0 2.5]);
grid on;
grid minor;
endFunction to Solve Conservative form of Equation
function [m_c, M_c, v_c, Temp_c, p_c, rho_c, Temp_throat_c, rho_throat_c, p_throat_c, M_throat_c] = conservative_form(x,dx,n,C,a,nt,gamma,throat)
%% Initialization step:
% Initializing density and temperature:
for i = 1:length(x)
if x(i) >= 0 && x(i) < 0.5
rho_c(i) = 1;
Temp_c(i) = 1;
else if x(i) >= 0.5 && x(i) < 1.5
rho_c(i) = 1.0 -0.366*(x(i) - 0.5);
Temp_c(i) = 1.0 - 0.167*(x(i) - 0.5);
else if x(i) >= 1.5 && x(i) <= 3.0
rho_c(i) = 0.634 - 0.3879*(x(i) - 1.5);
Temp_c(i) = 0.833 - 0.3507*(x(i) - 1.5);
end
end
end
end
% Initializing velocity and pressure:
v_c = 0.59./(rho_c.*a);
p_c = rho_c.*Temp_c;
% Initializing the solution vectors:
U1 = rho_c.*a;
U2 = rho_c.*a.*v_c;
U3 = rho_c.*a.*((Temp_c./(gamma - 1)) + (gamma/2)*v_c.^2);
% Time Iteration Loop
for k = 1:nt
% Preserving the old solution vectors:
U1_old = U1;
U2_old = U2;
U3_old = U3;
% Adaptive Time-step control:
dt = min(C.*dx./(sqrt(Temp_c) + v_c));
% Defining the flux vectors:
F1 = U2;
F2 = (U2.^2./U1) + ((gamma-1)/gamma)*(U3 - (gamma*U2.^2)./(2*U1));
F3 = (gamma*U2.*U3./U1) - (gamma*(gamma-1)/2)*(U2.^3./U1.^2);
%% Solution by Mc-Cormack Method:
for j = 2:n-1
% Defining the source term:
J2(j) = (1/gamma).*rho_c(j).*Temp_c(j).*((a(j+1)-a(j))/dx);
% Predictor step:
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 the new values of solution vectors:
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
% Updating the primitive variables and flux vectors:
rho_c = U1./a;
v_c = U2./U1;
Temp_c = (gamma -1)*(U3./U1 - (gamma/2)*v_c.^2);
p_c = rho_c.*Temp_c;
F1 = U2;
F2 = (U2.^2./U1) + ((gamma-1)/gamma)*(U3 - (gamma*U2.^2)./(2*U1));
F3 = (gamma*U2.*U3./U1) - (gamma*(gamma-1)/2)*(U2.^3./U1.^2);
for j = 2:n-1
% Updating the source term:
J2(j) = (1/gamma)*rho_c(j).*Temp_c(j)*((a(j)-a(j-1))/dx);
% Corrector step:
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
% Updating the solution vectors after corrector step:
% Averaging the predictor and corrector steps:
dU1dt = 0.5*(dU1dt_p + dU1dt_c);
dU2dt = 0.5*(dU2dt_p + dU2dt_c);
dU3dt = 0.5*(dU3dt_p + dU3dt_c);
% Calculating the new solution vectors:
for j = 2:n-1
U1(j) = U1_old(j) + dU1dt(j)*dt;
U2(j) = U2_old(j) + dU2dt(j)*dt;
U3(j) = U3_old(j) + dU3dt(j)*dt;
end
% Applying the boundary conditions:
% left side/ inlet conditions
U1(1) = rho_c(1).*Temp_c(1);
U2(1) = 2*U2(2) - U2(3);
U3(1) = U1(1)*(Temp_c(1)/(gamma-1)+(gamma/2)*v_c(1).^2);
% right side/outlet conditions:
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);
% Final updation of the primitive variables:
rho_c = U1./a;
v_c = U2./U1;
Temp_c = (gamma -1)*(U3./U1 - (gamma/2)*v_c.^2);
p_c = rho_c.*Temp_c;
%%
% Defining the Values of mach Number and mass flow rate:
M_c = v_c./sqrt(Temp_c);
m_c = rho_c.*a.*v_c;
% Calculating values of specific quantites at the throat:
M_throat_c(k) = M_c(throat); % Mach Number
p_throat_c(k) = p_c(throat); % Pressure Ratio
Temp_throat_c(k) = Temp_c(throat); % Non-dimensional Temperature
rho_throat_c(k) = rho_c(throat); % Density Ratio
% Comparison of non-dimensional mass-flow rates at different times:
mass2 = figure(4);
hold on;
if k == 100
plot(x,m_c,'b','linewidth',1.75);
else if k == 200
plot(x,m_c,'g','linewidth',1.75);
else if k == 700
plot(x,m_c,'r','linewidth',1.75);
end
end
end
end
title('Comparison of non-dimensional mass flow rates at different time-instants (Conservative Form)');
xlabel('x/L (Non-dimensional distance through the nozzle)');
ylabel('rhoVA/rho_{o}V_{o}A_{o}');
legend('100Deltat','200Deltat','700Deltat');
axis([0 3 0.4 1.6]);
grid on;
grid minor;
endObservations
1. Timewise variation of Pressure, Temperature, Velocity and Mach Number at Nozzle Throat
a. Non Conservative Governing Equation Solution
b. Conservative Governing Equation Solution
Explanation: There are distinctly 2 sections of the graph both in Non Conservative and conservative forms. The section of the graph before 300 seconds, where the solution oscillates rapidly and the section of the graph after 300 seconds where the solution stablises. Hence, it may be inferred that the initial oscillations in solution occur due to unsteady expansion and compression of waves propagating through the nozzle throat, in the transient state. However, after 300 seconds, the flow stablizes and the physical parameters reach steady state and becomes approximately equal to the analytical solution.
It must be taken into account as well, that the oscillations in the transient state are greater in case of non conservative solution, as compared to the conservative solution. The reason for this may be concluded as in non conservative equation, the solution is being done directly with the primitive variables- hence, at certain points in transient flow, the variations in primitive variables are quite significant, because of individual round off error of the primitve variables. However, in case of conservative flow, though the errors in the primitive variables vary significantly, the conservation variables remain relatively more stable. Hence, the solution in conservative form has less oscillations.
2. Variation of Mach No. Pressure, Density and Temperature along the length of the nozzle
a. Non Conservative Governing Equation Solution
b. Conservative Governing Equation Solution
Explanation: The initial observation in the plot is that, both the non conservative as well as the conservative governing equations manifest the same plot Then, the explanation proceeds to the nature of the individual plots. In case of Mach No. the plot follows the expected pattern as the Mach No. ranges from subsonic level at the entry of the nozzle to sonic level at the throat (at lenght=1.5, mach No. = 1) and increases to supersonic level at the exit of the diverging section of the nozzle.
A common nature is observed on the other hand for the pressure, temperature and density plots. They remain fairly constant as the flow converges, through the convergent section of the nozzle, and begins to drop as the flow starts expanding isentropically from the throat to the diverging section of the nozzle. The isentropic expansion of supersonic flow results in drop of static pressure and temperature along the flow field and hence, the pressure and temperature drops off along the length of the diverging section of the nozzle
3. Non Dimenisonal Mass Flowrate through the nozzle at Various Timesteps
a. Non Conservative Governing Equation Solution
b. Conservative Governing Equation Solution
Explanation: Both the graphs here manifest a common trend. The variations in mass flow rate at the initial timestep is pretty significant, owing to its transience, which gradually die down as the flow field reaches steady state. At 700 seconds, the flow stablizes and the non dimensional mass flow ratio reaches closes t to the analytical solution.
The standard explanation for the larger oscillations in solutions of non conservative form as compared to the conservative form holds good in this case as well. The round off error the primitve variables of Non conservative form is much more than the conservative variables used in solving conservative equations.
4. Comparison of non dimensional mass flowrates of Conservative and Non Conservative forms
Explanation: The above solution compares non dimensional mass flowrates at non-conservative, conservative and analytical modes of solutions. The point that needs to be noted here is that, much like the other non conservative modes of solution, the fluctuation is much greater in non conservative form than the conservative form. The reason for this phenomenon may be observed to be very explicit. In non conservative form, the solutions are done with respect to the primitive variables. Hence the round off errors of the individual primitive variables ,multiply against each other to produce a significant error in solution, that makes the solution plot oscillate more about the mean position. In conservative form, the conservative variables balance the error of the primitve variables as compared to each other, and hence the rounded off error is much less. Hence, the solution plot of the conservative equation is closer to the analytical solution plot, as compared to the non conservative solution.
5. Grid Independence Test
Explanation: The above image illustrates the values of non dimensional density, temperature and pressure at the throat of the nozzle. The purpose of doing the grid independence test is to manifest the fact, that the solution obtained for the physical parameters is independent of mesh size. This is particularly important because the same flow problem may yield results which are poles apart, inspite of having the same flow geometry, boundary conditions and initial conditions, because of the grid spacings. This will make the solution thermodynamically inconsistent. hence, it is of much importance , that the results obtained by solving the fluid flow equations are grid independent.
In the above image, it is observed that, when the grid size increases from 31 to 61, there is hardly any change in density, temperature or pressure, both in cases of conservative or non conservative form. Infact, in both the cases, it may be stated that the values of the parameters obtained from the non conservative as well as the conservative equations, are approximatelyt equal to the analytical solution.
6. Non Conservative vs Conservative- Minimum No. of simulation cycles
Explanation: The above image illustrates the variation of non dimensional mass flowrate at the throat of the nozzle, along the duration of 1400 seconds of flow. As expected, from the plot it can be concluded that the variation in non conservative solution is much more than conservative solution. The mathematical reason for this is that the partial derivatives of the primitive variables are multiplied against each other in the non conservative equations i.e. the round off errors of the primitve variables multiply against each other to enhance the magnitude of error, leading to higher oscillations in solution. Conservative equations are solved in terms of conservative variables, which do not multiply against each other and hence, the round off error created is minimal.
From the plot it is clearly visible that the conservative solution stabilizes after 520 time steps into a horizontal line, whereas the non conservative solution stabilizes only after 680 time steps. Hence, it may be inferred that the minimum number of time steps to convergence for conservative solution is 520 and the minimum number of time steps to convergence for non conservative solution is 680.
7. Non Conservative vs Conservative -The faster method
For 1400 time steps, with a nozzle legth of 3 units and a grid refinement of 31 grids, the computational time required to solve non conservative governing equation is 10.846 seconds while the time required to solve conservative governing equation is 7.7492 seconds. Hence, it may be concluded that the consevative mathod is the faster method to solution.