Simulation of a 1D Super-sonic nozzle flow simulation using McCormack Method
JORGE LOZANO
Introduction to CFD using Matlab and OpenFOAM
Simulation of a 1D Super-sonic nozzle flow simulation using McCormack Method
OBJECTIVE
The objective in this project is to write code solve the 1D supersonic nozzle flow equations using the McCormack Method:
- Implement both the conservative and non-conservative forms of the governing equations.
- Perform grid dependence test.
- With 1400 steps of time, and 31 nodes within the dimension of the nozzle, and 101 nodes to perform the grid dependency.
Steps to perform in this document:
- Write separate functions for conservative and non-conservative forms
- Figure out the minimum number of cycles for which the simulation should be run for convergence
- Compare the normalized mass flow rate between the conservative forms and non-conservative forms - make suitable conclusions
- Which method is faster?
Theory
In computational fluid dynamics, the McCormack method is a widely used discretization scheme for the numerical solution of hyperbolic partial differential equations. This second-order finite difference method was introduced by Robert W. McCormack in 1969.
The McCormack method is a variation of the two-step Lax–Wendorff scheme but is much simpler in application. To illustrate the algorithm, consider the following first order hyperbolic equation
The application of McCormack method to the above equation proceeds in two steps, a predictor step which is followed by a corrector step.
The McCormack method is well suited for nonlinear equations (Inviscid Burgers equation, Euler equations, etc.) The order of differencing can be reversed for the time step (i.e., forward/backward followed by backward/forward). For nonlinear equations, this procedure provides the best results. For linear equations, the McCormack scheme is equivalent to the Lax–Wendorff method.
Unlike first-order upwind scheme, the McCormack does not introduce diffusive errors in the solution. However, it is known to introduce dispersive errors (Gibbs phenomenon) in the region where the gradient is high.
The McCormack’s method is processed by applying an explicit scheme, of a predictor-corrector part each time step, according to the next two equations:
By using the governing equations, Euler’s 1D compressible equation applicable to a nozzle geometry, obtaining pressure, density, temperature and velocity respectively, are:
Mass equation:
Momentum equation:
Energy equation:
To obtain relative results of the variables according to the geometry x axis distance, the following equations are proposed to program:
- Non-dimensional nozzle length:
- Non-dimensional nozzle velocity (to obtain Mach’s number):
- Non-dimensional nozzle pressure:
- Non-dimensional nozzle temperature:
- Non-dimensional nozzle mass density:
Procedure
For the following case three analysis were made, by analytical solution, non-conservative data and conservative, a using functions with in MATLAB, a programing code was performed to obtain the results and data needed to conclude with the given data results, and the grid dependency effect within important variables.
Code for general 1D nozzle
close all
clear all
clc
%Non conservative form Validation data
%TIME STEP 1
RhoInitialnc=[1,0.955,0.927,0.900,0.872,0.844,0.817,0.789,0.760,0.731,0.701,0.670,0.637,0.603,0.567,0.531,0.494,0.459,0.425,0.392,0.361,0.330,0.301,0.271,0.242,0.213,0.184,0.154,0.125,0.095,0.066]
Tinitialnc=[1,0.972,0.950,0.929,0.908,0.886,0.865,0.843,0.822,0.800,0.778,0.755,0.731,0.707,0.682,0.656,0.631,0.605,0.581,0.556,0.533,0.510,0.487,0.465,0.443,0.421,0.398,0.376,0.354,0.332,0.309]
Vinitialnc=[0.111,0.212,0.312,0.411,0.508,0.603,0.695,0.784,0.870,0.954,1.035,1.113,1.188,1.260,1.328,1.394,1.455,1.514,1.568,1.619,1.666,1.709,1.748,1.782,1.813,1.838,1.858,1.874,1.884,1.890,1.895]
%TIME STEP 1400
Rhofinalnc=[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];
Tfinalnc=[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]
Vfinalnc=[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]
%Conservative form Validation data
%TIME STEP 1
RhoInitialc=[1,1,1,1,1,0.999,0.963,0.927,0.891,0.854,0.818,0.781,0.744,0.707,0.670,0.633,0.594,0.555,0.517,0.478,0.440,0.401,0.362,0.324,0.285,0.246,0.208,0.169,0.131,0.093,0.063]
Tinitialc=[1,1,1,1,1,1,0.983,0.967,0.950,0.934,0.917,0.900,0.883,0.866,0.849,0.833,0.800,0.766,0.731,0.695,0.660,0.625,0.590,0.554,0.519,0.484,0.448,0.412,0.375,0.324,0.200]
Vinitialc=[0.099,0.111,0.125,0.141,0.160,0.187,0.228,0.271,0.325,0.389,0.467,0.557,0.656,0.759,0.854,0.930,0.979,0.992,0.975,0.939,0.893,0.848,0.809,0.781,0.766,0.768,0.791,0.846,0.949,1.133,1.438]
%TIME STEP 1400
Rhofinalc=[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]
Tfinalc=[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]
Vfinalc=[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]
n=31;
m=0.579*ones(n)
x=linspace(0,3,n);
grid=linspace(0,3,31)
tic
[rho_nc,T_nc,p_nc,M_nc,m_nc,n_nc]= non_conservative_form;
figure(7)
plot(grid,RhoInitialnc,\'+\',grid,Tinitialnc,\'o\',grid,Vinitialnc,\'*\',grid,rho_nc,\'r\',grid,T_nc,\'b\',x,M_nc,\'g\')
xlabel(\'NOZZLE LENGTH\');
ylabel(\'MASS FLOW RATE\');
legend(\'RHO-INITIAL\',\'T-INITIAL\',\'V-INITIAL\',\'RHO-1\',\'T-1\',\'M-1\');
title(sprintf(\'FINAL VALUES RESULTS GRID DEPENDENCY NON-CONSERVATIVE\'));
figure(8)
plot(grid,Rhofinalnc,\'+\',grid,Tfinalnc,\'o\',grid,Vfinalnc,\'*\',grid,rho_nc,\'r\',grid,T_nc,\'b\',x,M_nc,\'g\')
xlabel(\'NOZZLE LENGTH\');
ylabel(\'NON-DIMENSIONAL VALUES\');
legend(\'RHO-FINAL\',\'T-FINAL\',\'V-FINAL\',\'RHO-1400\',\'T-1400\',\'M-1400\');
title(sprintf(\'FINAL VALUES RESULTS GRID DEPENDENCY NON-CONSERVATIVE\'));
tic
[rho_c,T_c,p_c,M_c,m_c,n_c]= conservative_form;
figure(9)
plot(grid,RhoInitialc,\'+\',grid,Tinitialc,\'o\',grid,Vinitialc,\'*\',grid,rho_c,\'r\',grid,T_c,\'b\',grid,M_c,\'g\')
xlabel(\'NOZZLE LENGTH\');
ylabel(\'MASS FLOW RATE\');
legend(\'RHO-INITIAL\',\'T-INITIAL\',\'V-INITIAL\',\'RHO-1\',\'T-1\',\'M-1\');
title(sprintf(\'FINAL VALUES RESULTS GRID DEPENDENCY CONSERVATIVE\'));
figure(10)
plot(grid,Rhofinalc,\'+\',grid,Tfinalc,\'o\',grid,Vfinalc,\'*\',grid,rho_c,\'r\',grid,T_c,\'b\',grid,M_c,\'g\')
xlabel(\'NOZZLE LENGTH\');
ylabel(\'NON-DIMENSIONAL VALUES\');
legend(\'RHO-FINAL\',\'T-FINAL\',\'V-FINAL\',\'RHO-1400\',\'T-1400\',\'M-1400\');
title(sprintf(\'FINAL VALUES RESULTS GRID DEPENDENCY CONSERVATIVE\'));
figure(11)
plot(x,m_nc,\'g\',x,m_c,\'b\',x,m,\'r\');
xlabel(\'LENGTH\');
ylabel(\'MASS FLOW RATE\');
title(sprintf(\'GRID= %d\\n\\nMASS FLOW ANALYSIS COMPARISON\',n));
legend(\'NON-CONSERVATION\',\'CONSERVATION\',\'ANALYTICAL\');Code for conservative solution
function [rho,t,p,M,m,count]= conservative_form
n=31;
nt=1400;
C=0.5;
tol=1e-4;
x=linspace(0,3,n);
dx=x(2)-x(1);
A=1+2.2*(x-1.5).^2;
throat=(n+1)/2;
nt_x=linspace(1,nt,nt);
gamma=1.4;
error=1;
count=1;
m_old=0;
m=0.579*ones(n);
u2_i=0.59;
for i=1:n
if ((x(i)>=0)&&(x(i)<=0.5))
rho(i)=1;
t(i)=1;
elseif ((x(i)<1.5)&&(x(i)>0.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
v(i)=u2_i/(rho(i)*A(i));
u1(i)=rho(i)*A(i);
u2(i)=rho(i)*A(i)*v(i);
u3(i)=rho(i)*A(i)*(t(i)/(gamma-1)+(gamma/2)*v(i)^2);
end
u1_initial=u1;
u2_initial=u2;
u3_initial=u3;
for k=1:n
dt_x(k)=C*dx/(t(k)^0.5 +v(k));
end
dt=min(dt_x);
for j=1:nt
u1_old=u1;
u2_old=u2;
u3_old=u3;
for i=2:n-1
f1_i=u2(i);
f1_ip1=u2(i+1);
f2_i=u2(i)^2 / (u1(i)) + ((gamma-1)/(gamma))*(u3(i)-(gamma/2)*u2(i)^2/u1(i));
f2_ip1=u2(i+1)^2 / (u1(i+1)) + ((gamma-1)/(gamma))*(u3(i+1)-(gamma/2)*u2(i+1)^2/u1(i+1));
j2=(1/gamma)*rho(i)*t(i)*((A(i+1)-A(i))/dx);
f3_i=gamma*u2(i)*u3(i)/u1(i) - (gamma*(gamma-1)/2)*u2(i)^3/u1(i)^2;
f3_ip1=gamma*u2(i+1)*u3(i+1)/u1(i+1) - (gamma*(gamma-1)/2)*u2(i+1)^3/u1(i+1)^2;
du1_dtp(i)=-(f1_ip1-f1_i)/dx;
du2_dtp(i)=j2-(f2_ip1-f2_i)/dx;
du3_dtp(i)=-(f3_ip1-f3_i)/dx;
u1(i)=u1(i)+du1_dtp(i)*dt;
u2(i)=u2(i)+du2_dtp(i)*dt;
u3(i)=u3(i)+du3_dtp(i)*dt;
end
rho=u1./A;
v=u2./u1;
t=(gamma-1)*((u3./u1)-(gamma/2)*v.^2);
for i=2:n-1
f1_i=u2(i);
f1_ic1=u2(i-1);
f2_i=u2(i)^2 / (u1(i)) + ((gamma-1)/(gamma))*(u3(i)-(gamma/2)*u2(i)^2/u1(i));
f2_ic1=u2(i-1)^2 / (u1(i-1)) + ((gamma-1)/(gamma))*(u3(i-1)-(gamma/2)*u2(i-1)^2/u1(i-1));
j2=(1/gamma)*rho(i)*t(i)*((A(i)-A(i-1))/dx);
f3_i=gamma*u2(i)*u3(i)/u1(i) - (gamma*(gamma-1)/2)*u2(i)^3/u1(i)^2;
f3_ic1=gamma*u2(i-1)*u3(i-1)/u1(i-1) - (gamma*(gamma-1)/2)*u2(i-1)^3/u1(i-1)^2;
du1_dtc(i)=-(f1_i-f1_ic1)/dx;
du2_dtc(i)=j2-(f2_i-f2_ic1)/dx;
du3_dtc(i)=-(f3_i-f3_ic1)/dx;
end
du1_dt=0.5*(du1_dtp+du1_dtc);
du2_dt=0.5*(du2_dtp+du2_dtc);
du3_dt=0.5*(du3_dtp+du3_dtc);
for i=2:n-1
u1(i)=u1_old(i)+du1_dt(i)*dt;
u2(i)=u2_old(i)+du2_dt(i)*dt;
u3(i)=u3_old(i)+du3_dt(i)*dt;
end
u1(1)=rho(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);
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);
for i=1:n
rho(i)=u1(i)/A(i)
v(i)=u2(i)/u1(i)
t(i)=(gamma-1)*(u3(i)/u1(i)-gamma/2*v(i)^2)
end
M=v./t.^0.5;
p=rho.*t;
m=rho.*v.*A;
m_t(j,:)=rho.*v.*A;
rho_throat(j)=rho(throat);
T_throat(j)=t(throat);
p_throat(j)=p(throat);
M_throat(j)=M(throat);
error=max(abs(m-m_old));
if error>tol
count=count+1;
m_old=m;
end
end
time_c=toc;
figure(4);
subplot(4,1,1);
plot(x,M,\'r\');
ylabel(\'Mach No.\');
title(sprintf(\'GRID= %d\\n\\nNOZZLE RESULTS CONSERVATION \\n COMPUTATIONAL TIME=%d \\n MACH NUMBER\',n,time_c));
subplot(4,1,2);
plot(x,rho,\'r\');
ylabel(\'p/po\');
title(sprintf(\'NON-DIMENSIONAL DENSITY\'));
subplot(4,1,3);
plot(x,t,\'r\');
ylabel(\'T/To\');
title(sprintf(\'NON-DIMENSIONAL TEMPERATURE\'));
subplot(4,1,4);
plot(x,p,\'r\');
xlabel(\'NOZZLE LENGTH\');
ylabel(\'P/Po\');
title(sprintf(\'NON-DIMENSIONAL PRESSURE\'));
figure(5)
plot(nt_x,rho_throat,\'g\');
hold on
plot(nt_x,T_throat,\'y\');
plot(nt_x,p_throat,\'r\');
plot(nt_x,M_throat,\'b\');
title(sprintf(\'GRID= %d\\n\\nFLOW FIELD VARIABLES: CONSERVATION NOZZLE\',n));
xlabel(\'TIME STEPS\');
legend(\'DENSITY\',\'TEMPERATURE\',\'PRESSURE\',\'MACH NUMBER\');
figure(6)
plot(x,u2_initial,\'--\',x,m_t(50,:),x,m_t(100,:),x,m_t(150,:),x,m_t(200,:),x,m_t(700,:),x,m,\':\')
xlabel(\'NOZZLE LENGTH: NON-DIMENSIONAL\');
ylabel(\'MASS FLOW: NON-DIMENSIONAL\');
legend(\'TIME STEP 1\',\'TIME STEP 50\',\'TIME STEP 100\',\'TIME STEP 150\',\'TIME STEP 200\',\'TIME STEP 700\',\'ANALYTICAL SOLUTION\');
title(sprintf(\'GRID= %d\\n\\nNON-DIMENSIONAL MASS FLOW: CONSERVATION NOZZLE \"m\" \',n));
endCode for non-conservative solution
function [rho,t,p,M,m,count]= non_conservative_form
n=31;
nt=1400;
x=linspace(0,3,n);
dx=x(2)-x(1);
A=1 + 2.2*(x-1.5).^2;
throat=(n+1)/2;
nt_x=linspace(1,nt,nt);
gamma=1.4;
error=1;
count=1;
m_old=0;
C=0.5;
tol=1e-4;
rho=1-0.3146*x;
t=1-0.2314*x;
v=(0.1+1.09.*x).*t.^(1/2);
m=0.579*ones(n);
m_initial=rho.*A.*v;
for k=1:n
dt_x(k)=C*dx/(t(k)^(1/2)+v(k));
end
dt=min(dt_x);
for j=1:nt
rho_old=rho;
v_old=v;
t_old=t;
for i=2:n-1
drho_dx=(rho(i+1)-rho(i))/dx;
dv_dx=(v(i+1)-(v(i)))/dx;
dt_dx=(t(i+1)-(t(i)))/dx;
dlogA_dx=(log(A(i+1))-(log(A(i))))/dx;
drho_dtp(i)=-(rho(i)*dv_dx)-(rho(i)*v(i)*dlogA_dx)-(v(i)*drho_dx);
dv_dtp(i)=-(v(i)*dv_dx)-(1/gamma)*(dt_dx+(t(i)/rho(i))*drho_dx);
dt_dtp(i)=-(v(i)*dt_dx)-(gamma-1)*t(i)*(dv_dx+v(i)*dlogA_dx);
rho(i)=rho_old(i)+drho_dtp(i)*dt;
v(i)=v_old(i)+dv_dtp(i)*dt;
t(i)=t_old(i)+dt_dtp(i)*dt;
end
for i=2: n-1
drho_dx=(rho(i)-rho(i-1))/dx;
dv_dx=(v(i)-v(i-1))/dx;
dt_dx=(t(i)-t(i-1))/dx;
dlogA_dx=(log(A(i))-log(A(i-1)))/dx;
drho_dtc(i)=-(rho(i)*dv_dx)-(rho(i)*v(i)*dlogA_dx)-(v(i)*drho_dx);
dv_dtc(i)=-(v(i)*dv_dx)-(1/gamma)*(dt_dx+(t(i)/rho(i))*drho_dx);
dt_dtc(i)=-(v(i)*dt_dx)-(gamma-1)*t(i)*(dv_dx+v(i)*dlogA_dx);
end
drho_dt=0.5*(drho_dtp+drho_dtc);
dv_dt=0.5*(dv_dtp+dv_dtc);
dt_dt=0.5*(dt_dtp+dt_dtc);
for i=2: n-1
rho(i)=rho_old(i)+drho_dt(i)*dt
v(i)=v_old(i)+dv_dt(i)*dt
t(i)=t_old(i)+dt_dt(i)*dt
end
v(1)=2*v(2)-v(3);
rho(n)=2*rho(n-1)-rho(n-2);
v(n)=2*v(n-1)-v(n-2);
t(n)=2*t(n-1)-t(n-2);
M=v./t.^0.5;
p=rho.*t;
m=rho.*v.*A;
m_t(j,:)=rho.*v.*A;
rho_throat(j)=rho(throat);
T_throat(j)=t(throat);
p_throat(j)=p(throat);
M_throat(j)=M(throat);
error=max(abs(m-m_old));
if error>tol
count=count+1;
m_old=m;
end
end
time_nc=toc;
figure(1);
subplot(4,1,1);
plot(x,M,\'g\');
ylabel(\'Mach No.\');
title(sprintf(\'GRID= %d\\n\\nNOZZLE RESULTS NON-CONSERVATION \\n COMPUTATIONAL TIME=%d \\n MACH NUMBER\',n,time_nc));
subplot(4,1,2);
plot(x,rho,\'g\');
ylabel(\'p/po\');
title(sprintf(\'NON-DIMENSIONAL DENSITY\'));
subplot(4,1,3);
plot(x,t,\'g\');
ylabel(\'T/To\');
title(sprintf(\'NON-DIMENSIONAL TEMPERATURE\'));
subplot(4,1,4);
plot(x,p,\'g\');
xlabel(\'NOZZLE LENGTH\');
ylabel(\'P/Po\');
title(sprintf(\'NON-DIMENSIONAL PRESSURE\'));
figure(2)
plot(nt_x,rho_throat,\'g\');
hold on
plot(nt_x,T_throat,\'y\');
plot(nt_x,p_throat,\'r\');
plot(nt_x,M_throat,\'b\');
title(sprintf(\'GRID= %d\\n\\nFLOW FIELD VARIABLES: NON-CONSERVATION NOZZLE \',n))
xlabel(\'TIME STEPS\');
legend(\'DENSITY\',\'TEMPERATURE\',\'PRESSURE\',\'MACH NUMBER\');
figure(3)
plot(x,m_initial,\'--\',x,m_t(50,:),x,m_t(100,:),x,m_t(150,:),x,m_t(200,:),x,m_t(700,:),x,m,\'r\');
xlabel(\'NOZZLE LENGTH: NON-DIMENSIONAL\');
ylabel(\'MASS FLOW: NON-DIMENSIONAL\');
legend(\'TIME STEP 1\',\'TIME STEP 50\',\'TIME STEP 100\',\'TIME STEP 150\',\'TIME STEP 200\',\'TIME STEP 700\',\'ANALYTICAL SOLUTION\');
title(sprintf(\'GRID= %d\\n\\nNON-DIMENSIONAL MASS FLOW: NON CONSERVATION NOZZLE \"m\" \',n));
endAccording to the Figure 1., by comparing the grid dependency within the system, it appears to perform very similar, taking notice that the courant number is maintained the same, as well as the other variables. These considerations are processed while a non-conservative scheme is analyzed.
Taking more computational time with more grid nodes. Having less computational time with 31 grid nodes, but more when it has 101 grid nodes.
Figure 1. Nozzle results for non conservation scheme; (a) grid=31, (b) grid=101.
In figure 2., we can conclude according to the flow filed variables vs. the time steps simulated that, the courant number and the corrections part of the method must be modified with more grid nodes, due to the convergance of these important results, seeing a convergance time step of 500 with 31 grid nodes, and 1400 with 101 nodes, in a non-conservative system.
Figure 2. Flow field variables for non conservation nozzle scheme; (a) grid=31, (b) grid=101.
In figure 3.,are obtained the mass-flow variations along the complete simulation with 31 and 101 grid nodes, concluding that with a non-conservation nozzle analyzis it generates more discrepancy in the beginning of the simulation with more nodes, but converges in both to obtaing a general valid range of results with a non-dimensional nozzle length.
Figure 3. Mass flow results for non conservation nozzle scheme; (a) grid=31, (b) grid=101.
According to the Figure 4., by comparing the grid dependency within the system, it appears to perform very similar, taking notice that the courant number is maintained the same, as well as the other variables. These considerations are processed while a conservative scheme is analyzed.
Taking more computational time with more grid nodes. Having more computational time with 31 grid nodes, but less when it has 101 grid nodes.
Figure 4. Flow field variables for conservation nozzle scheme; (a) grid=31, (b) grid=101.
In figure 5., we can conclude according to the flow filed variables vs. the time steps simulated that, the courant number and the corrections part of the method must be modified with more grid nodes, due to the convergance of these important results, seeing a convergance time step of 600 with 31 grid nodes, and 1400 with 101 nodes, in a conservative system.
Figure 5. Flow field variables for conservation nozzle scheme; (a) grid=31, (b) grid=101.
In figure 6.,are obtained the mass-flow variations along the complete simulation with 31 and 101 grid nodes, concluding that with a non-conservation nozzle analyzis it generates more discrepancy in the end of the simulation with more nodes, but has more convergence in both obtaing a general valid range of results with a dimensional nozzle length.
Figure 6. Mass flow results for conservation nozzle scheme; (a) grid=31, (b) grid=101.
The following Figure 7., shows a grid comparison within the different schemes used and the mass-flow in the results, having more variation in the non-conservation cheme while having more grid nodes, also, it obtains a more closed result with the conservation scheme relatively to the analytical solution which it maintains the same.
Figure 7. Mass flow results comparison of schemes analyzed; (a) grid=31, (b) grid=101.
With Figure 8., and Figure 9., the initial and final conditions and results are obtained, compared with a similar nozzle geometry, although by the Mach results within the four (initial and final, conservative and non-conservative schemes) that the geometry after the middle section is modified, obtating a different set of values, by simulating different geometrys and conditions the range of velocity obtained are necessary and importatn for engineering design with in the nozzle of air-vehicles affected by these results. Having good and expected results for density and tempereture conditions.
Figure 8. Final values grid dependency at initial conditions; (a) non-conservative, (b) conservative.
Figure 9. Final values grid dependency at final conditions; (a) non-conservative, (b) conservative.
Conclusion
A Supersonic nozzle 1D simulation was performed with different grid node numbers, obtaining a comparison and a relationship with 31 nodes, along the given data, was compared to obtain a different geometry effect of the nozzle seen in the Mach’s number obtained (M>1, supersonic system), all cases with a non-conservative, and a conservative scheme, obtaining the computational time taken and comparing the general equation applications with the John Anderson document.
The non-conservative scheme seems to perform faster in general although the conservative performs better while converging in relation to the analytical solution. A more accurate courant value could be proposed for more grid nodes in relation to the mass flow data comparisons obtained but variating very little with in the relative non-dimensional results, having more accuracy.