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

% Function code for conservative form (C)

function [mass_flow_rate_C, pressure_C, mach_number_C, rho_C, v_C, T_C, mass_flow_rate_throat_C, pressure_throat_C, mach_number_throat_C,rho_throat_C, v_throat_C,T_throat_C] = conservative_form(x,dx,n,gamma,nt,C)

% initial profile of density and tempreature

for i=1:n
    if x(i)>=0 && x(i)<0.5
        rho_C(i)=1;
        T_C(i)=1;
    elseif x(i)>=0.5 && x(i)<1.5
        rho_C(i)=1-0.366*(x(i)-0.5);
        T_C(i)=1-0.167*(x(i)-0.5);
    elseif x(i)>=1.5 && x(i)<=3
        rho_C(i)=0.634-0.3879*(x(i)-1.5);
        T_C(i)=0.833-0.3507*(x(i)-1.5);
    end
end

% profile
a = 1+2.2*(x-1.5).^2;
throat = find(a==1);
v_C = 0.59./(rho_C.*a);
pressure_C = (rho_C.*T_C);

U1 = rho_C.*a;
U2 = rho_C.*a.*v_C;
U3 = rho_C.*a.*(T_C./(gamma-1))+((gamma/2).*v_C.^2);

for k=1:nt
    dt = min(C.*dx./(sqrt(T_C)+v_C));
    %dt = 0.0207;
    U1_old = U1;
    U2_old = U2;
    U3_old = U3;

    % calculating flux vector
    F1 = U2;
    F2 = ((U2.^2)./U1)+((gamma-1)/gamma)*(U3-(gamma/2)*((U2.^2)./(U1)));
    F3 = ((gamma*U2.*U3)./U1) - (gamma*(gamma-1)*((U2.^3)./(2*U1.^2)));

    % predictor method
    for j = 2:n-1
        j2(j)=((1/gamma)*rho_C(j)*T_C(j)*((a(j+1)-a(j))/dx));
        % contiunity equation
        dU1dt_p(j) = -((F1(j+1)-F1(j))/dx);
        % momentum equation
        dU2dt_p(j) = -((F2(j+1)-F2(j))/dx)+j2(j);
        % energy equation
        dU3dt_p(j) = -((F3(j+1)-F3(j))/dx);
        
        % new values
        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

    % update flux vector
    F1 = U2;
    F2 = ((U2.^2)./U1)+((gamma-1)/gamma)*(U3-(((gamma/2)*U2.^2)./(U1)));
    F3 = ((gamma*U2.*U3)./U1) - (((gamma/2)*(gamma-1)*(U2.^3)./(U1.^2)));

    % corrector method
    for j=2:n-1
        j2(j)=((1/gamma)*rho_C(j)*T_C(j)*((a(j)-a(j-1))/dx));

        % equations 
        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

    % compute avg 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
    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

    % boundary condition
    % at inlet
    U1(1) = rho_C(1)*a(1);
    U2(1) = 2*U2(2)-U2(3);
    U3(1) = (U1(1)*((T_C(1)/(gamma-1))+((gamma/2)*(v_C(1)).^2)));
    v_C(1) = U2(1)./U1(1);

    % at 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);
    % final value update
    rho_C = U1./a;
    v_C = U2./U1;
    T_C = (gamma-1)*((U3./U1-(gamma/2)*(v_C).^2));

    mass_flow_rate_C = rho_C.*a.*v_C;
    pressure_C = rho_C.*T_C;
    mach_number_C = (v_C./sqrt(T_C));
    
    rho_throat_C(k) = rho_C(throat);
    T_throat_C(k) = T_C(throat);
    v_throat_C(k) = v_C(throat);
    mass_flow_rate_throat_C(k) = mass_flow_rate_C(throat);
    pressure_throat_C(k) = pressure_C(throat);
    mach_number_throat_C(k) = mach_number_C(throat);

    % plotting
    % comparison of non dimensional mass flow rate at different timestep
    figure(4)
    if k==50
        plot(x,mass_flow_rate_C, 'k', 'linewidth',1.2);
        hold on 
    elseif k==100
        plot(x,mass_flow_rate_C, 'g','linewidth',1.2);
        hold on
    elseif k==200
        plot(x,mass_flow_rate_C, 'm','linewidth',1.2);
        hold on
    elseif k==400
        plot(x,mass_flow_rate_C, 'y','linewidth',1.2);
        hold on
    elseif k==800
        plot(x,mass_flow_rate_C, 'c','lineWidth',1.2);
        hold on
    end

end
  % labelling 
    title('mass flow rate at different timestep (C)')
    xlabel('distance')
    ylabel('mass flow rate')
    legend({'50t^h timestep','100^t^h timestep','200^t^h timestep','400^t^h timestep','800^t^h timestep'})
    axis([0 3 0.4 1])
    grid on
end

% Function code for non conservative form (NC)

function [mass_flow_rate_NC, pressure_NC, mach_number_NC, rho_NC, v_NC, T_NC, mass_flow_rate_throat_NC, pressure_throat_NC, mach_number_throat_NC,rho_throat_NC, v_throat_NC,T_throat_NC] = non_conservative_form(x,dx,n,gamma,nt,C)



% Initial boundary condition
rho_NC = 1-0.3146*x;
T_NC = 1-0.2314*x;
v_NC = (0.1+1.09*x).*T_NC.^0.5;
a = 1+2.2*(x-1.5).^2;
throat = find(a==1);

dt = 0.0207;

for k = 1:nt
    dt = min(C.*dx./(sqrt(T_NC)+v_NC));
    rho_old = rho_NC;
    v_old = v_NC;
    T_old = T_NC;
    % predictor method
    for j=2:n-1
        drhodx = (rho_NC(j+1)-rho_NC(j))/dx;
        dvdx = (v_NC(j+1)-v_NC(j))/dx;
        dTdx = (T_NC(j+1)-T_NC(j))/dx;
        dlogadx = (log(a(j+1))-log(a(j)))/dx;
        % continuity equation'
        drhodt_p(j) = -rho_NC(j)*dvdx-rho_NC(j)*v_NC(j)*dlogadx-v_NC(j)*drhodx;
        % momentum equation
        dvdt_p(j) = -v_NC(j)*dvdx -(1/gamma)*(dTdx+(T_NC(j)*drhodx/rho_NC(j)));
        % energy equation
        dTdt_p(j)= -v_NC(j)*dTdx-(gamma-1)*T_NC(j)*(dvdx+v_NC(j)*dlogadx);

        %solution update
         rho_NC(j) = rho_NC(j)+drhodt_p(j)*dt;
         v_NC(j) = v_NC(j) + dvdt_p(j)*dt;
         T_NC(j) = T_NC(j) + dTdt_p(j)*dt;
    end
    % corrector method
    for j=2:n-1
    
        drhodx = (rho_NC(j)-rho_NC(j-1))/dx;
        dvdx = (v_NC(j)-v_NC(j-1))/dx;
        dTdx = (T_NC(j)-T_NC(j-1))/dx;
        dlogadx = (log(a(j))-log(a(j-1)))/dx;

        % continuity equation'
        drhodt_c(j) = -rho_NC(j)*dvdx-rho_NC(j)*v_NC(j)*dlogadx-v_NC(j)*drhodx;
        % momentum equation
        dvdt_c(j) =-v_NC(j)*dvdx -(1/gamma)*(dTdx+(T_NC(j)*drhodx/rho_NC(j)));
        % energy equation
        dTdt_c(j) = -v_NC(j)*dTdx-(gamma-1)*T_NC(j)*(dvdx+v_NC(j)*dlogadx);
    end
    % compute the 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
    for i =2:n-1
        rho_NC(i) = rho_old(i)+drhodt(i)*dt;
        v_NC(i) = v_old(i)+ dvdt(i)*dt;
        T_NC(i) = T_old(i)+ dTdt(i)*dt;
    end

   % boundary condition
   % inlet 
   rho_NC(1) =1;
   v_NC(1)= 2*v_NC(2)-v_NC(3);
   T_NC(1)= 1;

   % outlet
  
   rho_NC(n) = 2*rho_NC(n-1)-rho_NC(n-2);
   v_NC(n) = 2*v_NC(n-1)-v_NC(n-2);
   T_NC(n) = 2*T_NC(n-1)-T_NC(n-2);

   mass_flow_rate_NC = rho_NC.*a.*v_NC;
   pressure_NC = rho_NC.*T_NC;
   mach_number_NC = (v_NC./sqrt(T_NC));

   rho_throat_NC(k) = rho_NC(throat);
   T_throat_NC(k) = T_NC(throat);
   v_throat_NC(k) = v_NC(throat);
   mass_flow_rate_throat_NC(k) = mass_flow_rate_NC(throat);
   pressure_throat_NC(k) = pressure_NC(throat);
   mach_number_throat_NC(k) = mach_number_NC(throat);

   % plotting 
   % comparison of non dimensional mass flow rate at diffrent timestep
   figure(1)
   if k==50
       plot(x,mass_flow_rate_NC, 'r', 'linewidth',1.2);
       hold on
   elseif k==100
       plot(x, mass_flow_rate_NC,'g', 'linewidth', 1.2 );
       hold on
   elseif k==200
       plot(x, mass_flow_rate_NC,'b', 'linewidth', 1.2 );
       hold on
   elseif k==400
      plot(x, mass_flow_rate_NC,'c', 'linewidth', 1.2 );
      hold on
   elseif k==800
       plot(x, mass_flow_rate_NC,'m', 'linewidth', 1.2 );
       hold on
   end

end
   title('mass flow rate different timesteps (NC)')
   xlabel('distance')
   ylabel('mass flow rate')
   legend({'50t^h timestep','100^t^h timestep','200^t^h timestep','400^t^h timestep','800^t^h timestep'})
   axis([0 3 0 2])
   grid on
end

clear all
close all
clc

% input
n=31;
x=linspace(0,3,n);
dx = x(2)-x(1);
gamma = 1.4;
nt = 1400;
C=0.5;

% non conservative form
tic;
[mass_flow_rate_NC, pressure_NC, mach_number_NC, rho_NC, v_NC, T_NC, mass_flow_rate_throat_NC, pressure_throat_NC, mach_number_throat_NC,rho_throat_NC, v_throat_NC,T_throat_NC] = non_conservative_form(x,dx,n,gamma,nt,C)
time_elapsed_NC = toc;
fprintf('time elapsed for non conservative form: %0.4f', time_elapsed_NC);
hold on
% plotting
figure(2)
subplot(4,1,1)
plot(linspace(1,nt,nt), rho_throat_NC,'m')
ylabel('density')
legend({'density at throat'});
axis([0 1400 0 1.3]);
grid minor;
title('Timestep variation at throat (NC)')
subplot(4,1,2)
plot(linspace(1,nt,nt),T_throat_NC,'color','c')
ylabel('temperature')
legend({'temperature at throat'});
axis([0 1400 0.6 1]);
grid minor;
subplot(4,1,3)
plot(linspace(1,nt,nt),pressure_throat_NC,'color','g')
ylabel('pressure')
legend({'pressure at throat'});
axis([0 1400 0.4 1.1]);
grid minor;
subplot(4,1,4)
plot(linspace(1,nt,nt),mach_number_throat_NC,'color','y')
ylabel('mach number')
xlabel('timesteps')
legend({'mach number at throat'});
axis([0 1400 0.6 1.4]);
grid minor;

% steady state simulation for primitive variable for NC
figure(3);
subplot(4,1,1)
plot(x,rho_NC,'color','m')
ylabel('density')
legend({'density'});
axis([0 3 0 1]);
grid minor;
title('variation of flow rate distribution (NC)')
subplot(4,1,2)
plot(x,T_NC,'color','c')
ylabel('temperature')
legend({'temperature'});
axis([0 3 0 1]);
grid minor;
subplot(4,1,3)
plot(x,pressure_NC,'color','g')
ylabel('pressure')
legend({'pressure'});
axis([0 3 0 1]);
grid minor;
subplot(4,1,4)
plot(x,mach_number_NC,'color','y')
ylabel('mach number')
xlabel('nozzle length')
legend({'mach number'});
axis([0 3 0 4]);
grid minor;

% conservative form
tic;
[mass_flow_rate_C, pressure_C, mach_number_C, rho_C, v_C, T_C, mass_flow_rate_throat_C, pressure_throat_C, mach_number_throat_C,rho_throat_C, v_throat_C,T_throat_C] = conservative_form(x,dx,n,gamma,nt,C);
time_elapsed_C = toc;
fprintf('time elapsed for conservative form: %0.4f', time_elapsed_C);
% plotting 
figure(5)
subplot(4,1,1)
plot(rho_throat_C,'color','m')
ylabel('density')
legend({'density at throat'});
axis([0 1400 0 1]);
grid minor;
title('timestep variation at throat (C)')
subplot(4,1,2)
plot(T_throat_C,'color','C')
ylabel('temperature')
legend({'temperature at throat'});
axis([0 1400 0.6 1]);
grid minor;
subplot(4,1,3)
plot(pressure_throat_C,'color','g')
ylabel('pressure')
legend({'pressure at throat'});
axis([0 1400 0.2 0.8]);
grid minor;
subplot(4,1,4)
plot(mach_number_throat_C,'color','y')
ylabel('mach number')
xlabel('timestep')
legend({'mach number at throat'});
axis([0 1400 0.6 1.4]);
grid minor;

% steady state simulation for primitive variable for C
figure(6);
subplot(4,1,1)
plot(x,rho_C,'color','m')
ylabel('density')
legend({'density'});
axis([0 3 0 1]);
grid minor;
title(' variation of flow rate distibution (C)')
subplot(4,1,2)
plot(x,T_C,'color','c')
ylabel('temperature')
legend({'temperature'});
axis([0 3 0 1]);
grid minor;
subplot(4,1,3)
plot(x,pressure_C,'color','g')
ylabel('pressure')
legend({'pressure'});
axis([0 3 0 1]);
grid minor;
subplot(4,1,4)
plot(x,mach_number_C,'color','y')
ylabel('mach number')
xlabel('nozzle length')
legend({'mach number'});
axis([0 3 0 4]);
grid minor;

% plot for comparison of normalized mass flow rate

figure(7)

hold on
plot(x,mass_flow_rate_NC,'r','linewidth',1.2)
hold on
plot(x,mass_flow_rate_C,'g','linewidth',1.2)
hold on
grid on

legend('non conservative', 'conservative');
xlabel('length of domain')
ylabel('mass flow rate')
title('comparision of normalized mass flow rate')

The boundary condition are taken form the book Computational Fluid Dynamic by John D Anderson. We consider the steady isentropic flow through CD nozzle where the flow at inlet to the nozzle comes from a reservoir with cross sectional area of the reservoir is large (A tends to infinity) and velocity is very small (V tends to 0). The flow expand isentropically to supersonic speed at nozzle exit it is subsonic in the convergent, sonic at throat and supersonic in divergent section. Sonic flow at throat means that the velocity at this location is equal to the speed of sound.

Output: -

Output of the code has been uploaded as screen shot below. For variation of flow rate distribution graph we can observe that for both form they are giving same result for all the values. For time step variation of primitive at throat area we can observe that non conservative form for certain period take greater amount of time to achieve stable solution as comapred to conservative form.

Grid dependency is a factor to determine how problem can be solved faster so if we increase grid point we can observe that the values of variable doesnt change. We can conclude after observing all the result that conservative form converges slightly faster than non conservative form offering stable solution.