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);

%final Solution updates
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;
            dlogadx=(log(a(j+1))-log(a(j)))/dx;
            dadx=(a(j+1)-a(j))/dx;
            J2=(1/gamma).*rho_i(j).*t_i(j).*dadx;
            
            
            %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;
            dlogadx=(log(a(j))-log(a(j-1)))/dx;
            dadx=(a(j)-a(j-1))/dx;
            J2=(1/gamma).*rho_i(j).*t_i(j).*dadx;
            
            %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.

 

 

 

 

 

 


Projects by Tejas Shankhpal

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

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

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

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

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

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

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

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

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


Loading...

The End