1D Supersonic Nozzle flow simulation using Macormack method
Objective:
The main objective of this project is to write code solve the 1D supersonic nozzle flow equations using the Macormack Method. Before proceeding into this project first let us understand some basic things about nozzle.
Nozzle:
A nozzle is a device designed to control the direction or characteristics of a fluid flow (especially to increase velocity) as it exits (or enters) an enclosed chamber or pipe.
A nozzle is often a pipe or tube of varying cross sectional area, and it can be used to direct or modify the flow of a fluid (liquid or gas). Nozzles are frequently used to control the rate of flow, speed, direction, mass, shape, and/or the pressure of the stream that emerges from them. In a nozzle, the velocity of fluid increases at the expense of its pressure energy.
The pressure drops in a convergent nozzle because of the Bernoulli Principle. The Bernoulli principle states that an increase in the speed of a fluid occurs simultaneously with a decrease in the pressure exerted by the fluid.
Convergent Nozzle:
A convergent nozzle is a nozzle that starts big and gets smaller-a decrease in cross-sectional area. As a fluid enters the smaller cross-section, it has to speed up due to the conservation of mass. To maintain a constant amount of fluid moving through the restricted portion of the nozzle, the fluid must move faster.
The energy to make this fluid speed up has to come from somewhere. Some energy is in the random motion of molecules, which we observe as pressure. The energy in this random motion is converted into faster forward motion, known as stream flow. This change makes the pressure drop.
Pressure is inversely proportional to area, if everything else stays the same. In this case, the energy that causes pressure is converted to another type of energy, so both pressure and area decrease.
Here we consider the steady, isentropic flow through a convergent divergent nozzle as shown in the below fig,
The flow at the inlet to the nozzle comes from the reservoir where the pressure and the temperature are denoted by `P_0 and T_0` respectively. The cross sectional area of the reservoir is large (theorectically, `A-> oo`), and hence velocity is very small (`V->0`). Thus , `P_0 and T_0` are stagnation value ot total pressure and total temperature respectively. The flow expands isentropically to supersonic speeds at the nozzle exit, where the exit Pressure, Temperature, Velocity and Mach number are denoted by `P_e,T_e,V_e and M_e` respectively. The flow is locally subsonic in the convergent section of the nozzle, sonic at the throat (minimum area) and supersonic at the divergent section. The sonic flow (M=1) at the throat means that the local velocity at this location is equal to the local speed of sound. Using an asterisk to denote sonic flow values, we have at the throat V = V* = a*. Similarly the sonic flow values of the pressure and the temperature are denoted by p* and T* respectively. The area of the sonic throat are denoted by A*.
We assume that at a given section, where the cross sectional area is A, the flow properties are uniform across the section. Hence although the area of the nozzle changes as the function of distance along the nozzle x, and therefore in reality the flow field is two-dimensional(the flow varies in two dimensional xy space), we make the assumption that the flow properties vary only with x; this is tantamount to assuming uniform flow properties across any given cross section.Such flow is defined as Quasi one dimensional flow.
MacCormack's Technique:
The MacCormack method is an explicit finite-difference technique which is second order accurate in both space and time. This method is a popular technique as it provides satisfactory results for many fluid flow applications.MacCormack method is a variant of the Lax- Wendroff but is much simpler in its application because there is no need to evaluate the second order time derivatives with was the case for Lax-Wendroff method..MacCormack method is a two-step predictor- corrector sequence is used with forward differences on the predictor and with rearward differences on the corrector is a second order accurate method.
As we assume that the flow field at each grid point in time t is known, the flow field variables at each grid point at time `t+ Delta t` can be calculated using maccormack method in a predictor-corrector sequence.
In maccormack method,this can be obtained from,
`rho_(i,j)^(t+Deltat) = rho_(i,j)^t + (frac{del rho}{del t})_(av)* Deltat`
where, `(frac{delrho}{delt})(av)= 0.5* [(frac{delrho}{delt})(predicted)+ (frac{delrho}{delt})_(corrected)]`
Predicted step:
All the flow variables at time t are known, and we determine a predicted value of density, `(barrho)^(t+Deltat)` from the first two terms of Taylor series as follows,
`(barrho)^(t+Deltat) = rho_(i,j)^t + (frac{delrho}{delt})_(i,j)^t * Deltat`
here, `(frac{delrho}{delt})_(i,j)^t` can be obtained from the continuity equation and replacing the spatial derivatives on the right side with forward differences.
Corrected step:
In the corrected step, we first obtain a predicted value of time derivative at time `t+Deltat`,`(frac{bar(delrho)}{delt})(i,j)^(t+Deltat)` (also known as corrected time derivative), by substituting the predicted values of density (`(barrho)(i,j)^(t+Deltat)`),velocity (`(barv)_(i,j)^(t+Deltat)`) into the right side of the continuity equation, also replacing the spactial derivatives with rearward differences.
This sequence is repected at all grid points to obtain the density throught the flow field at time `t+Deltat`.To calculate velocity(v), temperature(t) and energy(e),the same technique is used and utilizing the momemtum and energy equations to obtain the average time derivatives via the predictor-corrector sequence, using forward differences on the predictor and rearward differences on the corrector.
Non- conservative form
Governing Equations:
Continuity equation:
`frac{delrho}{delt}= -rho*frac{delv}{delx}-rho*v*frac{del(ln A)}{dx}-v*frac{delrho}{delx}`
Momemtum equation:
`frac{delv}{delt}=-v*frac{delv}{delx}-frac{1}{gamma}*(frac{delT}{delx}+frac{T}{rho}*frac{delrho}{delx})`
Energy Equation:
`frac{delT}{delt}=-v*frac{delT}{delx}-(gamma-1)T(frac{delv}{delx}+v*frac{del(lnA)}{delx})`
For our solution of the quasi one dimensional nozzle flow of a calorically perfect gas, temperature(T) is a fundamental variable, therefore it is convinient to deal with it as the primary dependent variable in the energy equation.For a calorically perfect gas, energy, `e=c_v*T`.
Here, Density(`rho`), velocity (v), Temperature(T), nozzle distance(x), pressure (P), crossectional area (A), time(t) are expressed in terms of Non-dimensional variables.
Conservative Form:
Governing Equation:
Continuity Equation:
`frac{del(rhoA)}{delt}= -frac{del(rhoAv)}{delx}`
where, we assume
solution vector, `U_1= rhoA`
flux vector, `F_1= rhoAv`
Momemtum Equation:
`frac{del(rhoAv)}{delt}+ frac{del(rhoAv^2+(frac{1}{gamma}*PA))}{delx}= frac{1}{gamma}*P*frac{delA}{delx}`
where, we assume
solution vector, `U_2= rhoAv`
flux vector, `F_2= rhoAv^2+ frac{1}{gamma}PA`
source term, `J_2= frac{1}{gamma}Pfrac{delA}{delx}`
Energy Equation:
`frac{del(rho(frac{e}{gamma-1}+frac{gamma}{2}v^2)A)}{delt}+frac{del(rho(frac{e}{gamma-1}+frac{gamma}{2}v^2)vA+PAv)}{delx}=0`
where, we assume
solution vector, `U_3= rho(frac{e}{gamma-1}+frac{gamma}{2}v^2)A`
flux vector, `F_3=rho(frac{e}{gamma-1}+frac{gamma}{2}v^2)vA+PAv`
Therefore, above governing equations can be rewritten as
`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 of governing equations, the dependent variables( the variables for which we directly obtain numbers) are not the primitive variables.Therefore to obtain the primitive variables from above governing equations,
`rho= frac{U_1}{A}, v= frac{U_2}{U_1}, P= rhoT`
`T=e=(gamma-1)(frac{U_3}{U_1}-frac{gamma}{2}v^2)`
Pure Form of Flux terms: The flux terms are represented in terms of solution vectors instead to primitive variables(`rho,v,T`) in order to obtain a stable solution.Therefore,
`F_1= U_2`
`F_2= frac{U_2^2}{U_1}+ frac{gamma-1}{gamma}(U_3-frac{gamma}{2}frac{U_2^2}{U_1})`
`F_3= gammafrac{U_2U_3}{U_1}-frac{gamma(gamma-1)}{2}frac{U_2^3}{U_1^2}`
Courant- Friedrichs-Lowry(CFL) criterion:
`(Deltat)_i^t = Cfrac{Deltax}{a_i^t+v_i^t`
where, `v` is the local flow velocity at a point in the flow and `a` is the local speed of sound
For a stable solution, `C<=1` , we assume C=0.5.
we consider the minimum time step calculated at each grid point,i.e.,
`Deltat= min(Deltat_1^t,Deltat_2^t,....,Deltat_N^t)`
Matlab coding:
% 1D supersonic Nozzle flow
close all
clear all
clc
No_of_grid=[31 61]; % Grid matrix a=31 b=61
nt = 1400; % No of time steps
for i=1:length(No_of_grid)
if No_of_grid(i)==31
n = No_of_grid(i);
x = linspace(0,3,n);
dx = x(2)-x(1);
[rho_nc_a,t_nc_a,v_nc_a,dt_nc_a,M_nc_a,P_nc_a,count_nc_a,time_nc,mfr_nc_a,drhodt_a,dvdt_a,dtdt_a]=nonconservative_form(n,nt,x,dx);
[rho_c_a,T_c_a,v_c_a,dt_c_a,M_c_a,P_c_a,count_c_a,time_c,mfr_c_a]= conservative_form(n,nt,x,dx);
else
if No_of_grid(i)==61
m = No_of_grid(i);
y = linspace(0,3,m);
dy = y(2)-y(1);
[rho_nc_b,T_nc_b,v_nc_b,dt_nc_b,M_nc_b,P_nc_b,count_nc_b,time_nc,mfr_nc_b,drhodt_b,dvdt_b,dtdt_b]=nonconservative_form(m,nt,y,dy);
[rho_c_b,T_c_b,v_c_b,dt_c_b,M_c_b,P_c_b,count_c_b,time_c,mfr_c_b]= conservative_form(m,nt,y,dy);
end
end
end
% Plottings
figure(3)
subplot(4,1,1)
plot(x,M_nc_a,'r',y,M_nc_b,'b.');
hold on
xlabel('Non- dimensional length')
ylabel('Mach no')
title('Non Dimensionl Mach no variation of Non-conservative Form')
legend('31 grids','61 grids','location','west')
subplot(4,1,2)
plot(x,P_nc_a,'r',y,P_nc_b,'b.');
hold on
xlabel('Non- dimensional length')
ylabel('Pressure ratio')
title('Non dimensional pressure ratio variation of Non-conservative Form')
legend('31 grids','61 grids','location','west')
subplot(4,1,3)
plot(x,rho_nc_a,'r',y,rho_nc_b,'b.');
hold on
xlabel('Non- dimensional length')
ylabel('Density ratio')
title('Non dimensional Density ratio variation of Non-conservative Form')
legend('31 grids','61 grids','location','west')
subplot(4,1,4)
plot(x,t_nc_a,'r',y,T_nc_b,'b.');
hold on
xlabel('Non- dimensional length')
ylabel('Temp ratio')
title('Non dimensional Temp ratio variation of Non-conservative Form')
legend('31 grids','61 grids','location','west')
figure(4)
plot(x,mfr_c_a,'color','r','linewidth',2)
hold on
plot(x,mfr_nc_a,'color','b','linewidth',2)
plot([0 3],[0.579 0.579],'color','m','linewidth',2)
xlabel('non dimensional length')
ylabel('non dimensional mass flow rate')
title('Mass flow rate comparision of conservative and Non conservative form')
legend('conservative form','Non conservative form','Exact analytical solution')
figure(5)
subplot(4,1,1)
plot(x,M_c_a,'r',y,M_c_b,'b.');
hold on
xlabel('Non- dimensional length')
ylabel('Mach no')
title('Non Dimensional Mach No variation of conservative Form')
legend('31 grids','61 grids','location','west')
subplot(4,1,2)
plot(x,P_c_a,'r',y,P_c_b,'b.');
hold on
xlabel('Non- dimensional length')
ylabel('Pressure ratio')
title('Non dimensional pressure ratio variation of conservative Form')
legend('31 grids','61 grids','location','west')
subplot(4,1,3)
plot(x,rho_c_a,'r',y,rho_c_b,'b.');
hold on
xlabel('Non- dimensional length')
ylabel('Density ratio')
title('Non dimensional Density ratio variation of conservative Form')
legend('31 grids','61 grids','location','west')
subplot(4,1,4)
plot(x,T_c_a,'r',y,T_c_b,'b.');
hold on
xlabel('Non- dimensional length')
ylabel('Temp ratio')
title('Non dimensional Temp ratio variation of conservative Form')
legend('31 grids','61 grids','location','west')
Function Non conservative form:
function [rho,t,v,dt,M,P,count,time_nc,mfr,drhodt,dvdt,dtdt]= nonconservative_form(n,nt,x,dx)
% calculate initial profiles
rho = 1-0.3146*x; % density
t = 1-0.2314*x; % t = temp
v = (0.1+1.09*x).*t.^0.5; % volume
a = 1+2.2*(x-1.5).^2; % area
gamma = 1.4;
C = 0.5; % CFL
tol = 1e-5; % tolerance value
count = 0;
mfr = rho.*v.*a; % mass flow rate
for i=1:length(x)
if x(i)==1.5
th=i;
end
end
tic;
% Outer time loop
for k = 1:nt
rho_old = rho;
v_old = v;
t_old = t;
mfr_old = mfr;
% To find CFL condition
for j= 1:n
d_t(j)= C*(dx/(t(j)^0.5+v(j)));
end
dt(k)= min(abs(d_t));
% predictor method
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;
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)*(dtdx + (t(j)/rho(j))*drhodx);
% Energy equation
dtdt_p(j) = -v(j)*dtdx - (gamma-1)*t(j)*(dvdx + v(j)*dlogadx);
% Solution update
rho(j) = rho(j) + drhodt_p(j)*dt(k);
v(j) = v(j) + dvdt_p(j)*dt(k);
t(j) = t(j) + dtdt_p(j)*dt(k);
end
% Corrector method
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;
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)*(dtdx + (t(j)/rho(j))*drhodx);
% Energy equation
dtdt_c(j) = -v(j)*dtdx - (gamma-1)*t(j)*(dvdx + v(j)*dlogadx);
end
% Compute the average time derivative
for j = 2:n-1
drhodt(j) = 0.5*(drhodt_p(j) + drhodt_c(j));
dvdt(j) = 0.5*(dvdt_p(j) + dvdt_c(j));
dtdt(j) = 0.5*(dtdt_p(j) + dtdt_c(j));
end
% Final solution update
for i = 2:n-1
rho(i) = rho_old(i) + drhodt(i)*dt(k);
v(i) = v_old(i) + dvdt(i)*dt(k);
t(i) = t_old(i) + dtdt(i)*dt(k);
end
% Apply Boundary condition
% 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);
% Mass flow rate at diff time steps
mfr = rho.*v.*a;
if n==31
figure(1)
hold on
if k==1
plot(x,mfr,'linewidth',2)
else
if k==50
plot(x,mfr,'linewidth',2)
else
if k==100
plot(x,mfr,'linewidth',2)
else
if k==200
plot(x,mfr,'linewidth',2)
else
if k==500
plot(x,mfr,'linewidth',2)
else
if k==700
plot(x,mfr,'linewidth',2)
xlabel('Non dimensional distance through nozzle')
ylabel('Non dimensional mass flow rate')
title('Non dimensional mfr variation through nozzle at different no of timesteps')
legend('1 timestep','50 time steps','100 time steps','200 time steps','500 time steps','700 time steps')
end
end
end
end
end
end
end
% No of Timesteps to converge
error_mfr = max(mfr-mfr_old);
if error_mfr>=tol
count = count +1;
end
M= v./(t.^0.5); % Mach no
P= rho.*t; % Pressure
rho_th(k)=rho(th); % Density at throat of nozzle
t_th(k)= t(th); % Temp at throat of nozzle
M_th(k)= M(th); % Mach no at throat of nozzle
P_th(k)= rho(th)*t(th); % Pressure at throat of nozzle
time_nc = toc;
end
end
Function for Conservative form:
function [rho,T,v,dt,M,P,count,time_c,mfr]= conservative_form(n,nt,x,dx)
%initial condition
A=1+2.2.*(x-1.5).^2; % Area
gamma=1.4;
C=0.5; % CFL
Tol=1e-5; % Tolerance
count=0;
for i=1:length(x)
if x(i)==1.5
th=i;
end
end
tic;
for i= 1:n
if x(i)>=0 && x(i)<=0.5
rho(i)=1;
T(i)=1;
else
if 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);
else
if 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
end
end
v= 0.59./(rho.*A); % Volume
e=T;
mfr= rho.*v.*A;
% solution vectors
U1= rho.*A;
U2= rho.*A.*v;
U3= rho.*(e/(gamma-1)+(gamma/2).*v.^2).*A;
for i=1:nt
U1_old= U1;
U2_old= U2;
U3_old= U3;
rho_old=rho;
v_old=v;
T_old=T;
mfr_old=mfr;
% To find CFL condition
for j= 1:n
d_t(j)= C*(dx/(T(j)^0.5+v(j)));
end
dt(i)= min(abs(d_t));
% Flux vectors pure form
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);
for j= 2:n-1
% Source term
J2(j)= (1/gamma)*(rho(j)*T(j))*((A(j+1)-A(j))/dx);
end
% predicted step
for j= 2:n-1
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;
U1(j)=U1(j)+ dU1dt_p(j)*dt(i);
U2(j)=U2(j)+ dU2dt_p(j)*dt(i);
U3(j)=U3(j)+ dU3dt_p(j)*dt(i);
end
% Update predicted values
rho = U1./A;
T= (gamma-1).*((U3./U1)-(gamma/2).*(U2./U1).^2);
e=T;
v= U2./U1;
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);
for j=2:n-1
J2(j)= (1/gamma)*(rho(j)*T(j))*((A(j)-A(j-1))/dx);
end
% Corrected step
for j= 2:n-1
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
% Avg time derivative
dU1dt_av= 0.5.*(dU1dt_p+dU1dt_c);
dU2dt_av= 0.5.*(dU2dt_p+dU2dt_c);
dU3dt_av= 0.5.*(dU3dt_p+dU3dt_c);
for j=2:n-1
U1(j)= U1_old(j)+ dU1dt_av(j)*dt(i);
U2(j)= U2_old(j)+ dU2dt_av(j)*dt(i);
U3(j)= U3_old(j)+ dU3dt_av(j)*dt(i);
end
% Boundary condition
U2(1)= 2*U2(2)-U2(3);
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 update of primitive terms
rho = U1./A;
T= (gamma-1).*((U3./U1)-(gamma/2).*(U2./U1).^2);
e=T;
v= U2./U1;
% No of timesteps taken to converge
mfr= rho.*v.*A;
error_mfr=max(mfr-mfr_old);
if error_mfr>=Tol
count = count +1;
end
% Mass flow rate at diff time steps
if n==31
figure(2)
hold on
if i==1
plot(x,mfr,'linewidth',2)
else
if i==50
plot(x,mfr,'linewidth',2)
else
if i==100
plot(x,mfr,'linewidth',2)
else
if i==200
plot(x,mfr,'linewidth',2)
else
if i==500
plot(x,mfr,'linewidth',2)
else
if i==700
plot(x,mfr,'linewidth',2)
xlabel('Non dimensional distance through nozzle')
ylabel('Non dimensional Mass flow rate')
title('Non dimensional mfr variation through nozzle at different no of timesteps')
legend('1 timestep','50 time steps','100 time steps','200 time steps','500 time steps','700 time steps')
end
end
end
end
end
end
end
M= v./(T.^0.5); % Mach no
P= rho.*T; % Pressure
rho_th(i)=rho(th); % Density at throat of nozzle
T_th(i)= T(th); % Temp at throat of nozzle
M_th(i)= M(th); % Mach no at throat of nozzle
P_th(i)= rho(th)*T(th); % Pressure at throat of nozzle
time_c = toc;
end
endPlottings:
Non-conservative form:
Instantaneous distribution of the non dimensional mass flow as a function of distance through the nozzle at six different times during the time marching approach to the steady state.
Here in the below plotting , the blue color curve with the time step value 1 has the variation of density, volume and area which pertains the initial conditions.The distorted sine like variation of this curve is simply the product of assumed initial values for density and volume with the specified parabolic variation of the nozzle area ratio.
After 50 time steps the mass flow distribution through the nozzle has changed considerably and after 100 time steps it canges radically. Then however in 200 timesteps the mass flow distribution is beginning to settle down and after 700 time steps it is straight, horizontal line across the graph. This says that the mass flow has converged to a constant steady state value throughout the nozzle.
Conservative form:
Variation of mass flow distribution through the nozzle at different times during the time marching process, solution of the conservation form of the governing equations. The line which is in light blue color at 700 times steps indicates the excat analytical value.
Steady state Mass flow variation between Non-conservative and Conservative form:
From the above plot, let us compare the steady state variation of mass flow obtained with non conservative and conservative forms of the governing equation after 1400 time steps.
Here we can see that the steady state mass flow distribution predicted by the conservative form is much more satisfactory than the non conservative form.
Conservative form gives the distribution that is much closer to being a constant, whereas non conservative form has a sizebale variation with some oscillations at the both inlet and outlet boundaries.
The solution obatined from the conservative form is much closer to the excat analytical solution of 0.579.
Grid Independence Plot:
Grid Independence test is a serious consideration in the CFD. Here we have checked this with 31 grids and also by 61 grids. By observing the plots below there is not much difference in the solution for both the grids. The 31 grid is much closer to our analytial solution value, so we can conclude that by saying as the grid independent varibale. In most of the cases we need not the excat accurate solution. Increasing the grid number also increases our simulation time.
Non-conservative form:
Conservative Form:
Conclusion:
- The conservation form yields a better mass flow distribution. The conservation form simply does the better job for conserving the mass.
- There is no clear superiority of either forms in terms of accuracy of the results. The non conservative forms seems to produce slightly more accurate results for the primitive variables and the conservative form seems to produce slightly more accurate results for flux variables.The results in either case are certainly satisfactory.
- Comparing the amount of calculational effort to achieve a solution , conservative form requires more work than the non conservative form. This is due to the need to decode the primitive variables from the flux variables, such deoding process is not needed when we are solving non conservative form.
- Also comparing the the no of times steps to converge, for the same grid vlaue 31, the non conservative take less steps compared to conservative form. But in grid value 61 it is viceversa.