## Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method

This problem is solved for two forms of governing equations that are Conservative and Non conservative form for two types of grids that is 31 & 61:

Theory:

Maccormack Method is used to solve the Non-Conservative and Conservative form of governing equations which is as shown below:

1.Non-Conservative:

Code for Non-Conservative form is given as:

%1D Super-Sonic Nozzle Flow in Non-Conservative Form
close all
clc
%inputs
L=input('Enter the length of the domain: ');
n=input('Enter the number of grid points: ');
gamma=input('Enter the value of gamma: ');
%Specifying grid
x=linspace(0,L,n);
dx=x(2)-x(1);

%initial Conditions
rho=1-(0.3146*x);
T=1-(0.2314*x);
v=(0.1+(1.09*x)).*sqrt(T);
a=1+2.2*((x-1.5).^2);
n_throat= (n+1)/2;
c=0.5;
%time step
nt=1400;
dt=min(c.*(dx./((T.^0.5)+v)));
tic
%initial Condition
for k=1:nt
iter=0;
rho_old=rho;
T_old=T;
v_old=v;
massflow_old=rho_old.*v_old.*a;

%predictor step
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;
dlog_a_dx=(log(a(i+1))-log(a(i)))/dx;

%continuity equation
drho_dt_p(i)= -rho(i)*dv_dx-rho(i)*v(i)*dlog_a_dx-v(i)*drho_dx;

%energy equation
dT_dt_p(i)=-v(i)*dT_dx-(gamma-1)*T(i)*(dv_dx+v(i)*dlog_a_dx);

%momentum equation
dV_dt_p(i)=-v(i)*dv_dx-(dT_dx+(T(i)/rho(i))*drho_dx)/gamma;

%solution update
rho(i)=rho_old(i)+ (drho_dt_p(i)*dt);
T(i)=T_old(i)+ dT_dt_p(i)*dt;
v(i)=v_old(i)+dV_dt_p(i)*dt;

end

%Corrector step
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;
dlog_a_dx=(log(a(i))-log(a(i-1)))/dx;

%continuity equation
drho_dt_c(i)= -rho(i)*dv_dx-rho(i)*v(i)*dlog_a_dx-v(i)*drho_dx;

%energy equation
dT_dt_c(i)=-v(i)*dT_dx-(gamma-1)*T(i)*(dv_dx+v(i)*dlog_a_dx);

%momentum equation
dV_dt_c(i)=-v(i)*dv_dx-(dT_dx+(T(i)/rho(i))*drho_dx)/gamma;

end

drho_dt_avg=0.5*(drho_dt_c+ drho_dt_p);
dT_dt_avg=0.5*(dT_dt_c+dT_dt_p);
dV_dt_avg=0.5*(dV_dt_c+dV_dt_p);

for i= 2:n-1
rho(i)=rho_old(i)+drho_dt_avg(i)*dt;
T(i)=T_old(i)+dT_dt_avg(i)*dt;
v(i)=v_old(i)+dV_dt_avg(i)*dt;
end

%applying Boundary conditions;
%inlet Condition
v(1)=2*v(2)-v(3);

%outlet Condition
rho(n)=2*rho(n-1)-rho(n-2);
T(n)=2*T(n-1)-T(n-2);
v(n)=2*v(n-1)-v(n-2);

%mach number
M=v./sqrt(T);

%mass flow rate
massflow=rho.*v.*a;
p=rho.*T;

%throat values
rho_throat=rho(n_throat);
T_throat=T(n_throat);
p_throat=p(n_throat);
M_throat=M(n_throat);
v_throat=v(n_throat);

error=max(abs(massflow-massflow_old));
figure(1);
if k==1
plot(x,massflow,'--k')
hold on
elseif k==50
plot(x,massflow,'c')
hold on
elseif k==100
plot(x,massflow,'g')
hold on
elseif k==200
plot(x,massflow,'b')
hold on
elseif k==800
plot(x,massflow,'r')
hold on
elseif k==1400
plot(x,massflow,'y')
hold on
title('Mass flow Rates(Non Conservative)')
xlabel('X/L')
ylabel('\rhoVA/\rho_{o}v{o}a{o}')
legend('50\Deltat','100\Deltat','150\Deltat','200\Deltat','800\Deltat','1400\Deltat')
saveas(gcf,'massflow rate non conservative','jpeg')
iter=iter+1;
end
end
figure(2)
plot(x,rho)
title('Plot of Density Ratio vs Distance')
xlabel('X/L')
ylabel('\rho/\rho_{o}')
saveas(gcf,'Plot of Density vs Distance non conservative','jpeg')
hold on

figure(3)
plot(x,v,'g')
title('Plot of Velocity Ratio vs Distance')
xlabel('X/L')
ylabel('v/v{o}')
saveas(gcf,'Plot of Velocity vs Distance non conservative','jpeg')
hold on

figure(4)
plot(x,T,'r')
title('Plot of Temperature Ratio vs Distance')
xlabel('X/L')
ylabel('T/T{o}')
saveas(gcf,'Plot of Temperature vs Distance non conservative','jpeg')
hold on

figure(5)
plot(x,real(M),'y')
title('Plot of Mach number vs Distance')
xlabel('X/L')
ylabel('M')
saveas(gcf,'Plot of Mach number vs Distance non conservative','jpeg')
hold on

figure(6)
plot(x,p,'b')
title('Plot of Pressure Ratio vs Distance')
xlabel('X/L')
ylabel('P')
saveas(gcf,'Plot of Pressure Ratio vs Distance non conservative','jpeg')

simulation_time=toc;

sprintf('Computational Results from Non-Conservative Method are: \n Computation Time is = %d s\n Density Ratio From Results = %d \n Temperature Ratio= %d \n Velocity Ratio=%d \n Mach No=%d \n Pressure Ratio=%d',simulation_time,rho_throat,T_throat,v_throat,M_throat,p_throat)


For Grid Point:31

Inputs:

Plots

A.Plot of Mass Flow Rate vs Distance of the Nozzle:

B.Plot of Density Ratio vs Distance of the Nozzle:

C.Plot of Velocity Ratio vs Distance of the Nozzle:

D.Plot of  Temperature Ratio vs Distance of the Nozzle:

E.Plot of  Pressure Ratio vs Distance of the Nozzle:

F.Plot of  MACH NO vs Distance of the Nozzle:

Outputs:

For Grid = 61

2.Conservative:

Code for Conservative form is given as:

For Grid Point:31

%1D Super-Sonic Nozzle Flow in Conservative Form
close all
clc
%inputs
L=input('Enter the length of the domain: ');
n=input('Enter the number of grid points: ');
gamma=input('Enter the value of gamma: ');

%Specifying grid
x=linspace(0,L,n);
dx=x(2)-x(1);
nt=1400;

%initial Conditions
a=1+2.2*((x-1.5).^2);

for i=1:n
if (x(i)>=0 && x(i)<0.5)
rho_i(i)=1;
t_i(i)=1;
else
if (x(i)>=0.5 && x(i)<=1.5)
rho_i(i)=1-(0.366*(x(i)-0.5));
t_i(i)=1-(0.167*(x(i)-0.5));
else
if (x(i)>=1.5 && x(i)<=3.5)
rho_i(i)=0.634-(0.3879*(x(i)-1.5));
t_i(i)=0.833-(0.3507*(x(i)-1.5));
end
end
end

end

%velocity
v=(0.59)./((rho_i).*(a));
throat= (n+1)/2;
c=0.5;

%time step
for i=1:n
del_t(i)=min(c*(dx/((t_i(i)^0.5)+v(i))));
end
dt=min(del_t(i));

%variables
p=rho_i.*t_i;
U1=rho_i.*a;
U2=rho_i.*a.*v;
U3=rho_i.*((t_i./(gamma-1))+((gamma/2).*(v.^2))).*a;
tic
for k=1:nt

U1_old=U1;
U2_old=U2;
U3_old=U3;

%flux Variables
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));

%predictor step
for j=2:n-1
dF1dx=(F1(j+1)-F1(j))/dx;
dF2dx=(F2(j+1)-F2(j))/dx;
dF3dx=(F3(j+1)-F3(j))/dx;

%continuity equation
dU1dt_p(j)=-dF1dx;

%momentum equation
dU2dt_p(j)=-dF2dx+J2;

%energy equation
dU3dt_p(j)=-dF3dx;

%solution update
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
%corrector step
rho_i=U1./a;
v=U2./U1;
t_i=(gamma-1)*((U3./U1)-((gamma/2)*((U2./U1).^2)));

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
dF1dx=(F1(j)-F1(j-1))/dx;
dF2dx=(F2(j)-F2(j-1))/dx;
dF3dx=(F3(j)-F3(j-1))/dx;

%continuity equation
dU1dt_c(j)=-dF1dx;

%momentum equation
dU2dt_c(j)=-dF2dx+J2;

%energy equation
dU3dt_c(j)=-dF3dx;

end
%avg derivative
dU1dt=0.5*(dU1dt_c+dU1dt_p);
dU2dt=0.5*(dU2dt_c+dU2dt_p);
dU3dt=0.5*(dU3dt_c+dU3dt_p);

%final solution update
for i=2:n-1
U1(i)=U1_old(i)+dU1dt(i)*dt;
U2(i)=U2_old(i)+dU2dt(i)*dt;
U3(i)=U3_old(i)+dU3dt(i)*dt;
end
%applying boundary conditions
%inlet
U1(1)=rho_i(1)*a(1);
U2(2)=rho_i(2)*a(2)*v(2);
U3(1)=U1(1)*((t_i(1)/(gamma-1))+((gamma/2)*(v(1)^2)));

%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);
rho_i=U1./a;
v=U2./U1;
t_i=(gamma-1)*((U3./U1)-((gamma/2)*((U2./U1).^2)));

%calculating solutions
%mass flow rate
m=rho_i.*a.*v;
%pressure
p=rho_i.*t_i;
%Mach no
mach=v./sqrt(t_i);

%throat values
rho_t=rho_i(throat);
t_t=t_i(throat);
p_t=p(throat);
mach_t=mach(throat);
m_t=m(throat);
v_t=v(throat);
figure(1);
if k==1
plot(x,m,'--k')
hold on
elseif k==50
plot(x,m,'c')
hold on
elseif k==100
plot(x,m,'g')
hold on
elseif k==200
plot(x,m,'b')
hold on
elseif k==800
plot(x,m,'r')
hold on
elseif k==1400
plot(x,m,'y')
hold on
title('Mass flow Rates(Conservative)')
xlabel('X/L')
ylabel('\rhoVA/\rho_{o}v{o}a{o}')
legend('50\Deltat','100\Deltat','150\Deltat','200\Deltat','800\Deltat','1400\Deltat')
saveas(gcf,'massflow rate non conservative','jpeg')
end

end
figure(2)
plot(x,rho_i)
title('Plot of Density vs Distance')
xlabel('X/L')
ylabel('\rho/\rho_{o}')
saveas(gcf,'Plot of Density vs Distance conservative','jpeg')
hold on

figure(3)
plot(x,v,'g')
title('Plot of Velocity vs Distance')
xlabel('X/L')
ylabel('v/v{o}')
saveas(gcf,'Plot of Velocity vs Distance conservative','jpeg')
hold on

figure(4)
plot(x,t_i,'r')
title('Plot of Temperature vs Distance')
xlabel('X/L')
ylabel('T/T{o}')
saveas(gcf,'Plot of Temperature vs Distance  conservative','jpeg')
hold on

figure(5)
plot(x,mach,'y')
title('Plot of Mach number vs Distance')
xlabel('X/L')
ylabel('M')
saveas(gcf,'Plot of Mach number vs Distance conservative','jpeg')
hold on

figure(6)
plot(x,p,'b')
title('Plot of Pressure ratio vs Distance')
xlabel('X/L')
ylabel('P/P{o}')
saveas(gcf,'Plot of Pressure vs Distance conservative','jpeg')
simulation_time=toc;

sprintf('Computational Results from Conservative Method are: \n Computation Time is = %d s\n Density Ratio From Results = %d \n Temperature Ratio= %d \n Velocity Ratio=%d \n Mach No=%d \n Pressure Ratio=%d',simulation_time,rho_t,t_t,v_t,mach_t,p_t)


Inputs:

Plots

A.Plot of Mass Flow Rate vs Distance of the Nozzle:

OUTPUT:

FOR GRID POINT:61

CONCLUSION FROM THE RESULTS:

• From the Results it can be seen that the values coming from the computational results are matching with the exact results in the Anderson Book.
• Grid Independence check is performed in the Both Conservative and Non Conservative part of the domain and results are not varying.
• Moreover the fluctuations in mass flow rate plot for conservative part is lesser than the non conservative part. Hence Conservative Part gives much satisfying results comparing to the non conservative part.
• The accuracy of the solution increases as we are increasing Grid Points.
• The rate of decay of residuals is much faster in non Conservative form.
• For same number of grid points Non Conservative solution yeilds better results compare to conservative results.

### Compact Notation Derivation for a simple Mechanism Tejas Shankhpal · 2020-03-03 19:08:46

Reaction Mechanism: CO+O2 rarr_(k_(r1))^(k_(f1)) CO2+O O + H2O rarr_(k_(r2))^(k_(f2)) OH + OH CO + OH rarr_(k_(r3))^(k_(f3)) CO2 + H H + O2 rarr_(k_(r4))^(k_(f4)) OH + O kf1 and kr1 represent the forward and backward rate coefficients respectively. The rat Read more

### Combustion Efficiency Calculation after Preheating Tejas Shankhpal · 2020-02-17 12:57:33

Objective: 1.What is the effect of preheating on the adiabatic flame temperature? Change pre-heating from 298 to 600K Explain your results. For this problem, Methane is considered as fuel and combusted with air. CH4+2(O2+3.76N2)=CO2+2H2O+7.52N2 Code for the pr Read more

### Transient simulation of flow over a throttle body Tejas Shankhpal · 2020-02-15 20:24:58

Geometery :Geometery is imported to Converge Cfd as stl input file. Simulation Time Parameters: Start Time:0 s End Time:0.01s Boundary Conditions: Inlet Pressure:150000Pa Outlet pressure: 100000Pa Elbow and Throttle Plate: Wall Type Boundary Condition(law of wa Read more

### Steady state simulation of flow over a throttle body Tejas Shankhpal · 2020-02-15 17:24:00

Geometery :Geometery is imported to Converge Cfd as stl input file. Simulation Time Parameters: Start Time:0 Cycle End Time:15000 cycles Boundary Conditions: Inlet Pressure:150000Pa Outlet pressure: 100000Pa Elbow and Throttle Plate: Wall Type Boundary Conditio Read more

### combustion in fluent Tejas Shankhpal · 2019-11-25 20:23:21

Possible types of Combustion in fluent:Computational Fluid dyanmics modelling of combustion calls upon the proper selection and implementation of model suitable to faithfully represent the complex physical and chemical phenomenon associated with any combustion process. Read more

### cyclone seperator challenge Tejas Shankhpal · 2019-11-25 19:17:52

Empirical models used for cyclone separator efficiency: 1.IOZIA AND LEITH MODEL: Iozia and Leith (1990) logistic model is a modified version of Barth (1956) model which is developed based on force balance. The model assumes that a particle carried by thevortex endur Read more

### Rayleigh taylor instability Tejas Shankhpal · 2019-11-25 17:12:25

1.Practical CFD models that have been based on the mathematical analysis of Rayleigh Taylor waves. Kelvin–Helmholtz instability: The theory predicts the onset of instability and transition to turbulent flow in fluids of different densities moving at various spe Read more

### Meshing of Cylinder using ANSA Tejas Shankhpal · 2019-11-19 10:37:37

1.Geomtery: Geometery is imported in the ANSA from open option from file menu. 2.Meshing Select PID option from the bottom part of the menu. Now Go to presentation parametres and check the checkmark on the Draw cons,mid mode, like in ENT mode. Click A Read more

### Gate Valve Parametric Study Tejas Shankhpal · 2019-11-18 15:00:55

1.Gate Valve: Gate Valve is a valve that opens up by lifting a barrier or gate out of the path of the fluid.Gate Valve requires very little space along the pipe axis and hardly restricts the fluid flow when the gate is fully opened. 2.Geometry & Material Propertie Read more

### CHT analysis on an exhaust manifold Tejas Shankhpal · 2019-11-12 21:12:47

What is Conjugate Heat Transfer? The term conjugate heat transfer is used to describe processes which involve variations of temperature within solids and fluids, due to thermal interaction between the solids and fluids. A typical example is the heating or cooling of Read more