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

SIMULATION OF 1D SUPERSONIC FLOW NOZZLE USING MACROMAC METHOD:

*NON CONSERVATION FORM OF GOVERNING EQUATION:

CONTINUITY EQUATION

MOMENTUM EQUATION

ENERGY EQUATION

All the variables are non dimensional form.

initial condition

CONSERVATION FORM OF GOVERNING EQUATION

COUNTIUNITY EQUATION

MOMENTUM EQUATION

ENERGY EQUATION

The above equation can be simplfy into the following equation

where u1,u2,u3 are called the dolution vector, and F1,F2,F3 are the flux vector and j2 is called the source term.

INITIAL CONDITION

MATLAB CODE

clear all

close all
clc
n=31;
x=linspace(0,3,n);
tol=1e-4;
m_analytical=0.579*ones(1,n);
[m_nc]=supersonic_non_conservative(n,tol);
[m_c]=supersonic_conservative(n,tol);
figure(5)
plot(x,m_nc,'linewidth',2)
hold on
plot(x,m_c,'linewidth',2)
hold on
plot(x,m_analytical,'linewidth',2)
title('mass flow rate vs x')
xlabel('x')
ylabel('mass flow rate')
legend('non conservative','conservative','analytical')

function for non conservative

function [m_mc]=supersonic_non_conservative(n,tol)
tic
%grid 
x=linspace(0,3,n);

dx=x(2)-x(1);
gamma=1.4;
%specific throat
throat =(n+1)/2;
m_analytical=0.579*ones(1,n);

c=0.5;
%initial profiles
rho=1-0.3146*x;
t=1-0.2314*x;  % temperature
v=(0.1+1.09*x).*t.^0.5;
rho_old=rho'
v_old=v;
t_old=t;
%area
a=1+2.2*(x-1.5).^2;
%time steps
nt=1400;
nt_step=linspace(1,nt,nt);
for p=1:n
    dt_a(p)=c*dx/(t(p)^0.5+v(p));
end
dt=min(dt_a);
%outer time loop
iter=1;
error_rho=9e9;
error_v=9e9;
error_t=9e9;
for k=1:nt
    rho_old=rho;
    v_old=v;
    t_old=t;

    %predictor method
    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;

        % countinuity equation
        drhodt_p(j)=-rho(j)*dvdx-rho(j)*v(j)*dlogadx-v(j)*drhodx;

        %momentumequation
        dvdt_p(j)=-v(j)*dvdx-(1/gamma)*(dtdx+(t(j)/rho(j))*drhodx);

        %energy equation
        dtdt_p(j)=-v(j)*dtdx-(gamma-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 method
     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;

        % countinuity equation
        drhodt_c(j)=-rho(j)*dvdx-rho(j)*v(j)*dlogadx-v(j)*drhodx;

        %momentumequation
        dvdt_c(j)=-v(j)*dvdx-(1/gamma)*(dtdx+(t(j)/rho(j))*drhodx);

        %energy equation
        dtdt_c(j)=-v(j)*dtdx-(gamma-1)*t(j)*(dvdx+v(j)*dlogadx);
     end

     % compute avg time deravative
     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
     m=v./t.^0.5;
     m_nc=rho.*a.*v;
     m_nct(k,:)=rho.*a.*v;
     error_rho=max(abs(rho-rho_old));
     error_v=max(abs(v-v_old));
     error_t=max(abs(t-t_old));

     if (error_rho>tol && error_v>tol && error_t>tol)
         iter=iter+1;
     end

     %apply boundary conditions
     %inlet
     v(1)=2*v(2)-v(3);
     rho(n)=2*rho(n-1)-rho(n-2);
     t(n)=2*t(n-1)-t(n-2);

end
rho_throat=rho(throat);
v_throat=v(throat);
m_throat=m(throat);
figure(1)
subplot(4,1,2)
plot(x,rho,'linewidth',2)

xlabel('nozzle length')
ylabel('density ratio')
subplot(4,1,3)
plot(x,v,'r','LineWidth',2)

xlabel('nozzle length')
ylabel('velocity ratio')
subplot(4,1,4)
plot(x,t,'g','LineWidth',2)

xlabel('nozzle length')
ylabel('temperature ratio')
subplot(4,1,1)
plot(x,m,'m','LineWidth',2)
title(sprintf('no of timestep=%d','no of iteration to steady=%d',k,iter))
xlabel('nozzle length')
ylabel('mach number')
subtitle('nonconservative form')

figure(3)
plot(x,m_nct(50,:),x,m_nct(100,:),x,m_nct(200,:),x,m_nct(300,:),x,m_nct(500,:),'o',x,m_analytical,'LineWidth',2)
xlabel('non dimensionless distance x')
ylabel('non dimensionless mass flow rate')
title('variation of mass flow rate (non conservative)')
legend('50dt','100dt','200dt','300dt','700dt','analytical')
time_nc=toc;
fprintf('time of simulation for nonconservative=%gsn',time_nc);
end




     

function for conservative

function [m_c]=supersonic_conservative(n,tol)
tic
%grid
x=linspace(0,3,n);
dx=x(2)-x(1);
gamma=1.4;
throat=(n+1)/2;
m_analytical=0.579*ones(1,n);
%initial profiles
for q=1:n
    if((x(q)>0) && (x(q)<0.5))
        rho(q)=1;
        t(q)=1;
    elseif ((x(q)<1.5) && (x(q)>0.5))
        rho(q)=1-0.366*(x(q)-0.5);
        t(q)=1-0.167*(x(q)-0.5);
    elseif ((x(q)>=1.5)&& (x(q)<=3.5))
        rho(q)=0.634-0.3879*(x(q)-1.5);
        t(q)=0.833-0.03507*(x(q)-1.5);
    end
end
%area
a=1+2.2*(x-1.5).^2;
v=0.59./(rho.*a);
u1=rho.*a;
u2=rho.*a.*v;
u3=rho.*(t./(gamma-1)+(gamma/2)*v.^2).*a;
%time step
nt=1400;
nt_dtep=linspace(1,nt,nt);
c=0.5;
for p=1:n
    dt_a(p)=c*dx/(t(p)^0.5+v(p));
end
dt=min(dt_a);
%outer time loop
iter=1;
error_rho=9e9;
error_v=9e9;
error_t=9e9;

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=(gamma-1)*((u3_old./u1_old)-(gamma/2)*v_old.^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);

    %predictor method
    for j=2:n-1
        j2(j)=(1/gamma)*rho(j)*t(j)*(a(j+1)-a(j))/dx;

        %continuity equation
        du1dt_p(j)=-(f1(j+1)-f1(j))/dx;

        %momentumequation
        du2dt_p(j)=j2(j)-(f2(j+1)-f2(j))/dx;

        %energy equation
        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=(gamma-1)*((u3./u1)-(gamma/2)*v.^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)));

    %corrector method
    for j=2:n-1
        j2(j)=(1/gamma)*rho(j)*t(j)*(a(j)-a(j-1))/dx;
        %continuty equation
        du1dt_c(j)=-(f1(j)-f1(j-1))/dx;

        %momentum equation
        du2dt_c(j)=j2(j)-(f2(j)-f2(j-1))/dx;

        %energy equation
        du3dt_c(j)=-(f3(j)-f3(j-1))/dx;

    end

    %compute avg time deravative
    du1dt=0.5*(du1dt_p+du1dt_c);
    du2dt=0.5*(du2dt_p+du2dt_c);
    du2dt=0.5*(du2dt_p+du2dt_c);

    %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
    %apply boundary condition
    %inlet
    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);
    %outlet
    u2(n)=2*u2(n-1)-u2(n-2);
    u1(n)=2*u1(n-1)-u1(n-2);
    u3(n)=2*u3(n-1)-u3(n-2);

    %updating values
    rho=u1./a;
    v=u2./u1;
    t=(gamma-1)*((u3./u1)-(gamma/2)*v.^2);
    m=v./t.^0.5;
    m_c=rho.*a.*v;
    m_nc(k,:)=rho.*a.*v;

    error_rho=max(abs(rho-rho_old));
    error_v=max(bas(v-v_old));
    error_t=max(abs(t-t_old));

    if (error_rho>tol && error_v>tol && error_t>tol)
        iter=iter+1;
    end
end
rho_throat=rho(throat);
v_throat=v(throat);
m_throat=m(throat);

figure(2)
subplot(4,1,2)
plot(x,rho,'linewidth',2)

xlabel('nozzle length')
ylabel('density ratio')
subplot(4,1,3)
plot(x,v,'r','LineWidth',2)

xlabel('nozzle length')
ylabel('velocity ratio')
subplot(4,1,4)
plot(x,t,'g','LineWidth',2)

xlabel('nozzle length')
ylabel('temperature ratio')
subplot(4,1,1)
plot(x,m,'m','LineWidth',2)
title(sprintf('no of timestep=%d','no of iteration to steady state=%d',k,iter))
xlabel('nozzle length')
ylabel('mach number')
subtitle('conservative form')

figure(4)
plot(x,m_nc(50,:),x,m_nc(100,:),x,m_nc(200,:),x,m_nc(300,:),x,m_nc(700,:),'o',x,m_analytical,'linewidth',2);
xlabel('non dimensional distance x')
ylabel('non dimensionless mass floe rate')
title('variation of mass flow rate (conservative)')
legend('50dt','100dt','200dt','300dt','700dt','analytical')
time_c=toc;
fprintf('time of simulation for conservative=%gsn',time_c);
end

PLOTS

PROFILE FOR FLOW VARIABLE FOR NON CONSERVATIVE FORM

PROFILE FOR FLOW VARIABLE FOR CONSERVATIVE FORM

Variable for mass flow rate at different time step

for non conservation

for conservation

convergence:

through the above two graph the converence time step can be analysed we can say that at 700 time step the mass flow rate correspond to the analytical solution that means it reach to steady state solution.

Comparision of mass flow rate between conservative and non conservative

from this we can conclude that mass flow rate in conservative show better result correspond to analytical solution.