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

Matlab Code for Conservative Form:

%Conservative form
clear all;
clc all;
close all;

n=61;
x=linspace(0,3,n);
dx=x(2)-x(1);
gaamma=1.4;
c=0.5; %cfl number
e=10e9;
tol=10e-4; %convergence criterion
a=1+2.2*(x-1.5).^2;

%Initial conditions

for i=1:n

    if (x(i)<=0.5)
    rho(i)=1;t(i)=1;
  
  elseif (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);
  
  elseif(x(i)>=1.5 && x(i)<=3)
    rho(i)=0.634-0.3879*(x(i)-1.5);
    t(i)=0.833-0.3507*(x(i)-1.5);
  
  end

end

v=0.59./(rho.*a);

U1=rho.*a;
U2=rho.*a.*v;
U3=rho.*(t*(1/(gaamma-1))+(0.5*gaamma*v.^2)).*a;

%time steps
nt=1400;

%loop to calculate from courant number
for z=1:n
dt_x(z)=c*dx/(t(z)^0.5+v(z));
dt=min(dt_x);
end
iter=0;
tic;

%outer time loop
%while e>tol
for k=1:nt 

    U1_old=U1;
    U2_old=U2;
    U3_old=U3;
    rho_old=U1_old./a;
    v_old=U2_old./U1_old;
    t_old=(gaamma-1)*((U3_old./U1_old)-(0.5*gaamma*(v_old.^2)));
    F1=U2;
    F2=((U2.^2)./U1)+(((gaamma-1)/gaamma).*(U3-((0.5*gaamma)*(U2.^2)./U1)));
    %F2=(rho.*a.*v.^2)+(1/gaamma)*rho.*a.*t;
    F3=((gaamma*U2.*U3)./U1)-((gaamma*(gaamma-1)*0.5)*((U2.^3)./(U1.^2)));
    %F3=(rho.*a.*v).*((t/(gaamma-1))+(0.5*gaamma*v.^2))+rho.*a.*v;
    
    %predictor step
   for j=2:n-1
      J2(j)=(1/gaamma)*rho(j)*t(j)*((a(j+1)-a(j))/dx);
      %continuity eqn
      dU1dt_p(j)=-(F1(j+1)-F1(j))/dx;
      %momentum eqn
      dU2dt_p(j)=J2(j)-(F2(j+1)-F2(j))/dx;
      %energy eqn
      dU3dt_p(j)=-(F3(j+1)-F3(j))/dx;
      
      %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
   
    rho=U1./a;
    v=U2./U1;
    t=(gaamma-1)*((U3./U1)-(0.5*gaamma*(v.^2)));
    F1=U2;
    F2=((U2.^2)./U1)+(((gaamma-1)/gaamma).*(U3-((0.5*gaamma)*(U2.^2)./U1)));
    %F2=(rho.*a.*v.^2)+(1/gaamma)*rho.*a.*t;
    F3=((gaamma*U2.*U3)./U1)-((gaamma*(gaamma-1)*0.5)*((U2.^3)./(U1.^2)));
    %F3=(rho.*a.*v).*((t/(gaamma-1))+(0.5*gaamma*v.^2))+rho.*a.*v;
    
   
    %corrector step
    
    for i=2:n-1
      J2(i)=(1/gaamma)*rho(i)*t(i)*((a(i)-a(i-1))/dx);
      %continuity eqn
      dU1dt_c(i)=-(F1(i)-F1(i-1))/dx;
      %momentum eqn
      dU2dt_c(i)=J2(i)-(F2(i)-F2(i-1))/dx;
      %energy eqn
      dU3dt_c(i)=-(F3(i)-F3(i-1))/dx;
      
      %solution update
      U1(i)=U1(i)+dU1dt_c(i)*dt;
      U2(i)=U2(i)+dU2dt_c(i)*dt;
      U3(i)=U3(i)+dU3dt_c(i)*dt;
   end
   
   %calculate average time derivative
   dU1dt=0.5*(dU1dt_p+dU1dt_c);
   dU2dt=0.5*(dU2dt_p+dU2dt_c);
   dU3dt=0.5*(dU3dt_p+dU3dt_c);
   
   %final solution update
    for k=2:n-1
    U1(k)=U1_old(k)+dU1dt(k)*dt;
    U2(k)=U2_old(k)+dU2dt(k)*dt;
    U3(k)=U3_old(k)+dU3dt(k)*dt;    
    end
    
    %apply boundary conditions
    %inlet
    rho(1)=U1(1)./a(1);
    v(1)=U2(1)./(rho(1).*a(1));
    U2(1)=2*U2(2)-U2(3);
    U3(1)=U1(1).*(t(1)*(1/(gaamma-1))+(0.5*gaamma*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);
    
    %updating values
    rho=U1./a;
    v=U2./U1;
    t=(gaamma-1)*((U3./U1)-(0.5*gaamma*(v.^2)));
    
    
    e_rho=max(abs(rho-rho_old));
    e_v=max(abs(v-v_old));
    e_t=max(abs(t-t_old));
    error_matrix=[e_rho e_v e_t];
    e=max(error_matrix);
    iter=iter+1;
end
toc;
simulation_time=toc;
    %plot(x,v);
    %title(sprintf(\'Time step=%d\',k))
    %pause(0.03)
    %mass flow rate
   m=(rho.*v.*a);
   plot(x,m);

Matlab Code Non-Conservative Form:

%Non-Conservative form
clear all;
clc all;
close all;

n=61;
x=linspace(0,3,n);
dx=x(2)-x(1);
gaamma=1.4;
c=0.5; %cfl number
e=10e9;
tol=10e-4; %convergence criterion

%Initial conditions
rho=1-0.3146.*x;
t=1-0.2314.*x;
v=(0.1+1.09.*x).*t.^0.5;
a=1+2.2*(x-1.5).^2;


%time steps
nt=1400;

%loop to calcualate from courant number
for z=1:n
dt_x(z)=c*dx/(t(z)^0.5+v(z));
dt=min(dt_x);
end
iter=0;
tic;
%outer time loop
%while e>tol
for k=1:nt 
    
    rho_old=rho;
    v_old=v;
    t_old=t;
    
    %predictor step
    
   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 eqn
      drhodt_p(j)=(-rho(j)*dvdx)-(rho(j)*v(j)*dlogadx)-(v(j)*drhodx);
    
      %momentum eqn
      dvdt_p(j)=-v(j)*dvdx-(1/gaamma)*(dtdx+(t(j)/rho(j))*drhodx);
    
      %energy eqn
      dtdt_p(j)=(-v(j)*dtdx)-(gaamma-1)*t(j)*(dvdx+v(j)*dlogadx);
    
      %solution update
      rho(j)=rho(j)+drhodt_p(j)*dt;
      v(j)=v(j)+dvdt_p(j)*dt;
      t(j)=t(j)+dtdt_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;
      dtdx=(t(j)-t(j-1))/dx;
      dlogadx=(log(a(j))-log(a(j-1)))/dx;
    
      %continuity eqn
      drhodt_c(j)=-rho(j)*dvdx-rho(j)*v(j)*dlogadx-v(j)*drhodx;
    
      %momentum eqn
      dvdt_c(j)=-v(j)*dvdx-(1/gaamma)*(dtdx+(t(j)/rho(j))*drhodx);
    
      %energy eqn
      dtdt_c(j)=(-v(j)*dtdx)-(gaamma-1)*t(j)*(dvdx+v(j)*dlogadx);
    
   end
      
   %calculate average time derivative
   drhodt=0.5*(drhodt_p+drhodt_c);
   dvdt=0.5*(dvdt_p+dvdt_c);
   dtdt=0.5*(dtdt_p+dtdt_c);
   
   %final solution update
    for i=2:n-1
    rho(i)=rho_old(i)+drhodt(i)*dt;
    v(i)=v_old(i)+dvdt(i)*dt;
    t(i)=t_old(i)+dtdt(i)*dt;
    end
    
    %apply 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);  
    t(n)=2*t(n-1)-t(n-2);
    e_rho=max(abs(rho-rho_old));
    e_v=max(abs(v-v_old));
    e_t=max(abs(t-t_old));
    error_matrix=[e_rho e_v e_t];
    e=max(error_matrix);
    iter=iter+1;
    
    %plot(x,v);
    %title(sprintf(\'Time step=%d\',k))
    %pause(0.03)
    
%end
end
toc;
simulation_time=toc;
 %mass flow rate
      m=(rho.*v.*a);
      plot(x,m);

Normalized Mass Flow Rates:

Conservative Form:

Non-Conservative Form:

Conclusion:

From the above results, it can be concluded that the results corresponding to conservative form is near to the analytical results and shows good agreement with the actual results

Number of iterations for Convergence:

Conservative form: For a residual of 10e-4, the number of iterations required to reach convergence is 182

Non-Conservative form:For a residual of 10e-4, number of iterations required to reach convergence is 276

Simulation Time:

Conservative form:12.074 sec

Non-Conservative form:14.6266 sec

This shows that for 1400 timesteps Solution by Conservative form is faster than Non-Conservative form

Grid Dependency Test:

Conservative form:

Results @ throat for a grid size of 31:

rho_throat:0.65238 temp_throat:0.84107 vel_throat:0.89692

Results @ throat for a grid size of 61:

rho_throat:0.63899 temp_throat:0.83566 vel_throat:0.91216

Non-Conservative form:

Results @ throat for a grid size of 31:

rho_throat:0.6386 temp_throat:0.83645 vel_throat:0.91402

Results @ throat for a grid size of 61:

rho_throat:0.63766 temp_throat:0.91382 vel_throat:0.83544

From the results, it can be seen that there is not much difference in the results for 31 grid points and 61 grid points for both conservative & non-conservative form. This confirms the grid dependency test i.e. in order to save computation time and cost, a grid with 31 nodes will do the same work as that of a grid with 61 nodes.