Project

Quasi one-dimensional analysis of subsonic-supersonic isentropic nozzle flow using MacCormack s method

GGaurav Katiyar19 Dec 2019

Numerical solution of quasi one-dimensional subsonic-supersonic isentropic nozzle flow

(I) Non-conservative form of governing equations

1. Grid independence study: It is necessary that the numerical solution of a physical problem remains unaffected by the number of grid points adopted for the purpose of numerical discretization. Hence performing grid independence study is an essential step in CFD.

For the purpose of establishing the grid independence of our solution as well as validating our numerical codes we will compare the numerical values (obtained from the numerical code) of the parameters, at the throat section, with their exact analytical values.

(a) N=31 grid points

$\"\"$

The above picture shows the values of non-dimensional parameters, namely: density, temperature, pressure and Mach number, at the throat section for 31 grid points.

The absolute error in non-dimensional density is=0.74%

The absolute error in non-dimensional temperature is=0.42%

The absolute error in non-dimensional pressure is=1.18%

The absolute error in Mach number is=0.063%

Note that the numerical values obtained for 31 grid points are themselves within 2% of the exact analytical values.

(b) N=61 grid points

$\"\"$

The above picture shows the values of non-dimensional parameters, namely: density, temperature, pressure and Mach number, at the throat section for 61 grid points.

The absolute error in non-dimensional density is=0.56%

The absolute error in non-dimensional temperature is=0.29%

The absolute error in non-dimensional pressure is=0.87%

The absolute error in Mach number is=0.03%

Hence we can see that even upon doubling the number of grid points, there is practically no change in the values of the non-dimensional variables. In other words, 31 grid points provide a good enough spatial resolution for the current problem. Henceforth all the results discussed for the non-conservative form of the governing equations, will be with N=31 grid points.

2. Plot of residuals vs time steps

$\"\"$

3. Variation of non-dimensional variables at the throat with time steps

$\"\"$

The above plot shows the variation of non-dimensional density, non-dimensional velocity and non-dimensional temperature at the nozzle throat with time steps. As the number of time steps increases beyond 500 there is negligible change in the variables. In other words it can be said that in terms of individual variables, the steady state is achieved after 500 time steps. However the residuals only drop by approximately three orders of magnitude over 500 time steps.

4. Variation of normalized mass flow rate with time steps

$\"\"$

The above plot shows the normalized mass flow rates throughout the domain at different time instants. It can be seen that the spatial variations in the mass flow rate die down with an increase in time steps and after 700 time steps, the mass fow rate achieves almost a constant value.

5. Steady state results

$\"\"$

The above figure displays the steady state profiles for non-dimensional density (a), non-dimensional velocity (b), non-dimensional temperature (c) and non-dimensional pressure (d). It can be inferred from the plots that along the length of the nozzle, the density decreases, velocity increases, temperature decreases and pressure decreases.

(II) Conservative form of governing equations

(a) N=31 grid points

$\"\"$

The above picture shows the values of non-dimensional parameters, namely: density, temperature, pressure and Mach number, at the throat section for 31 grid points.

The absolute error in non-dimensional density is=1.19%

The absolute error in non-dimensional temperature is=0.29%

The absolute error in non-dimensional pressure is=1.5%

The absolute error in Mach number is=0.16%

Note that the numerical values obtained for 31 grid points are themselves within 2% of the exact analytical values.

(b) N=61 grid points

$\"\"$

The above picture shows the values of non-dimensional parameters, namely: density, temperature, pressure and Mach number, at the throat section for 61 grid points.

The absolute error in non-dimensional density is=0.42%

The absolute error in non-dimensional temperature is=0.16%

The absolute error in non-dimensional pressureis=0.60%

The absolute error in Mach number is=0.23%

Hence we can see that even upon doubling the number of grid points, there is practically no change in the values of the non-dimensional variables. In other words, 31 grid points provide a good enough spatial resolution for the current problem. Henceforth all the results discussed for the conservative form of the governing equations, will be with N=31 grid points.

2. Plot of residuals vs time steps

$\"\"$

The above graph shows the variation of residuals with time steps. It can be seen that the residuals keep dropping with an increase in time steps. This is expected as the solution approaches steady state. By 1200 time steps, the residuals have dropped to approximately 10^-4, which is negligible. Hence in terms of the average time derivatives, steady state is achieved after 1200 time steps but the simulation is still run until 1400 time steps, just to be conservative. Additionally, upon comparing the above plot with the residual plot for non-conservative form of governing equations, it can be inferred that the rate of decay of residuals is faster for the non-conservative form. Also the steady state values of residuals for the non-conservative form are orders of magnitude smaller than the steady state values of residuals for the conservative form.

3. Variation of non-dimensional variables at the throat with time steps

$\"\"$

4. Variation of normalized mass flow rate with time steps

$\"\"$

5. Steady state results

$\"\"$

(III) Comparison between non-conservative form and conservative form of governing equations

1. Steady state normalized mass flow rate

$\"\"$

The above plot shows the steady state distribution of normalized mass flow rate for both non-conservative form and conservative form of governing equations and compares them with the exact analytical solution. The plot for non-conservative form exhibits oscillations near the entry and exit of the nozzle. The conservative form leads to a more accurate steady state variation of mass flow rate. Furthermore the steady state solution achieved by the conservative form is much closer to the exact analytical solution. In other words, the conservative form conserves mass flow rate more effectively as compared to the non-conservative form. This is because mass flow rate is a primary dependent variable in conservative form as opposed to the non-conservative form where it is a derived variable.

2. Speed

$\"\"$

According to the above result, the conservative form takes slightly longer to achieve steady state. This can be attributed to the extra efforts required, in case of conservative form, for the conversion of dependent variables to primitive variables and calculation of flux variables at each time step.

(IV) Numerical Code

1. Main program

clc;
clear all;
close all;
% Inputs
N=31;
L=3;
x=linspace(0,L,N);
dx=x(2)-x(1);
gamma=1.4;
C=0.5;
% Initial profiles for Temperature, density and area
rho=1-0.3146*x;
T=1-0.2314*x;
A=1+2.2*(x-1.5).^2;
% Non-conservative form
tic
[rho1,v1,T1,p1,dmdt1]=non_conservative_form(N,L,x,dx,gamma,C,rho,T,A);
time1=toc;
% Conservative form
tic
[rho2,v2,T2,p2,dmdt2]=conservative_form(N,L,x,dx,gamma,C,rho,T,A);
time2=toc;
% Calculation of steady state normalized mass flow rate
m1=rho1.*A.*v1;
m2=rho2.*A.*v2;
m_exact=0.579*ones(1,N);
% Output
fprintf(\'Time taken to achieve steady state solution using the non-conservative form=%f \\n\',time1);
fprintf(\'Time taken to achieve steady state solution using the conservative form=%f \\n\',time2);
% Plotting
% 1. SS normalized mass flow rate
figure(5)
plot(x,m1,\'-k\')
hold on;
plot(x,m2,\'--k\');
plot(x,m_exact,\'-*k\');
xlabel(\'x/L\');
ylabel(\'Normalized mass flow rate (\\rhoAv/\\rho_0A*a_0)\');
legend(\'non-conservative form\',\'conservative form\',\'exact analytical solution\');
% 2. SS profiles
% (i) Non-conservative form
figure(6)
subplot(2,2,1)
plot(x,rho1);
xlabel(\'x/L\');
ylabel(\'\\rho/\\rho_0\');
title(\'(a)\');
subplot(2,2,2)
plot(x,v1);
xlabel(\'x/L\');
ylabel(\'v/a_0\');
title(\'(b)\');
subplot(2,2,3)
plot(x,T1);
xlabel(\'x/L\');
ylabel(\'T/T_0\');
title(\'(c)\');
subplot(2,2,4)
plot(x,p1);
xlabel(\'x/L\');
ylabel(\'p/p_0\');
title(\'(d)\');
% (ii) Conservative form
figure(7)
subplot(2,2,1)
plot(x,rho2);
xlabel(\'x/L\');
ylabel(\'\\rho/\\rho_0\');
title(\'(a)\');
subplot(2,2,2)
plot(x,v2);
xlabel(\'x/L\');
ylabel(\'v/a_0\');
title(\'(b)\');
subplot(2,2,3)
plot(x,T2);
xlabel(\'x/L\');
ylabel(\'T/T_0\');
title(\'(c)\');
subplot(2,2,4)
plot(x,p2);
xlabel(\'x/L\');
ylabel(\'p/p_0\');
title(\'(d)\');
% 3. Normalized mass flow rate at different time steps
% (i). Non-conservative form
figure(8)
plot(x,dmdt1(:,2),\'-k\',x,dmdt1(:,3),\'-*k\',x,dmdt1(:,4),\'-ok\',x,dmdt1(:,5),\'-sk\',x,dmdt1(:,6),\'-dk\');
xlabel(\'x/L\');
ylabel(\'Normalized mass flow rate (\\rhoAv/\\rho_0A*a_0)\')
legend(\'50\\Deltat\',\'100\\Deltat\',\'150\\Deltat\',\'200\\Deltat\',\'700\\Deltat\');
% (ii). Conservative form
figure(9)
plot(x,dmdt2(:,2),\'-k\',x,dmdt2(:,3),\'-*k\',x,dmdt2(:,4),\'-ok\',x,dmdt2(:,5),\'-sk\',x,dmdt2(:,6),\'-dk\');
xlabel(\'x/L\');
ylabel(\'Normalized mass flow rate (\\rhoAv/\\rho_0A*a_0)\')
legend(\'50\\Deltat\',\'100\\Deltat\',\'150\\Deltat\',\'200\\Deltat\',\'700\\Deltat\');

2. Function for non-conservative form

function [rho,v,T,p,dmdt]=non_conservative_form(N,L,x,dx,gamma,C,rho,T,A)
% Initial conditions
v=(0.1+1.09*x).*(T.^0.5);
% Initial mass flow rate
ct=1;
dmdt(:,1)=rho.*A.*v;

% Time steps
nt=1400;

% Outer time loop
for k=1:1:nt
    % Calculation of time step using CFL criterion
    for i=2:1:N-1
        time_step(i-1)=C*dx/(T(i)^0.5+v(i));
    end
    dt=min(time_step);
    rho_old=rho;
    v_old=v;
    T_old=T;
    % Predictor step
    for j=2:1: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 equation
        drhodt_p(j)=-rho(j)*dvdx-rho(j)*v(j)*dlogAdx-v(j)*drhodx;
        % Momentum equation
        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 step
    for j=2:1: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 equation
        drhodt_c(j)=-rho(j)*dvdx-rho(j)*v(j)*dlogAdx-v(j)*drhodx;
        % Momentum equation
        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
    % Calculation of the average time derivatives
    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 j=2:1:N-1
        rho(j)=rho_old(j)+drhodt(j)*dt;
        v(j)=v_old(j)+dvdt(j)*dt;
        T(j)=T_old(j)+dTdt(j)*dt;
    end
    % Boundary conditions
    % Inlet
    v(1)=2*v(2)-v(3);
    % Outlet
    rho(N)=2*rho(N-1)-rho(N-2);
    v(N)=2*v(N-1)-v(N-2);
    T(N)=2*T(N-1)-T(N-2);
    % Calculation of pressure
    p=rho.*T;
    % Calculation of residuals
    v_residual(k)=max(abs(v-v_old)/dt);
    rho_residual(k)=max(abs(rho-rho_old)/dt);
    T_residual(k)=max(abs(T-T_old)/dt);
    % Calculation of variables at the throat section
    rho_throat(k)=rho(0.5*(N+1));
    v_throat(k)=v(0.5*(N+1));
    T_throat(k)=T(0.5*(N+1));
    % Calculation of mass flow rate
    if(k==50||k==100||k==150||k==200||k==700)
        ct=ct+1;
        dmdt(:,ct)=rho.*A.*v;
    end
    % Plotting
    % 1. Residuals vs time steps
    figure(1)
    semilogy(rho_residual,\'-k\');
    hold on;
    semilogy(v_residual,\'--k\');
    semilogy(T_residual,\'-.k\');
    xlabel(\'Number of time steps\');
    ylabel(\'Residuals\');
    legend(\'|d\\rho/dt|_{av}\',\'|dv/dt|_{av}\',\'|dT/dt|_{av}\');
    % 2. Variables vs time steps
    figure(2)
    plot(rho_throat,\'-k\');
    hold on;
    plot(v_throat,\'--k\');
    plot(T_throat,\'-.k\');
    hold off;
    xlabel(\'Number of time steps\');
    ylabel(\'Non-dimensional variables at the throat\');
    legend(\'\\rho/\\rho_0\',\'v/a_0\',\'T/T_0\');
end
% Output
fprintf(\'Non-conservative form \\n \\n\');
fprintf(\'Non-dimensional density at the throat is=%f \\n\',rho(0.5*(N+1)));
fprintf(\'Exact analytical value of non-dimensional density is=0.634 \\n\');
fprintf(\'Non-dimensional temperature at the throat is=%f \\n\',T(0.5*(N+1)));
fprintf(\'Exact analytical value of non-dimensional temperature is=0.833 \\n\')
fprintf(\'Non-dimensional pressure at the throat is=%f \\n\',p(0.5*(N+1)));
fprintf(\'Exact analytical value of non-dimensional pressure is=0.528 \\n\');
fprintf(\'Mach number at the throat is=%f \\n\',v(0.5*(N+1))/T(0.5*(N+1))^0.5);
fprintf(\'Exact analytical value of Mach number is=1 \\n \\n\');
end

3. Function for conservative form

function [rho,v,T,p,dmdt]=conservative_form(N,L,x,dx,gamma,C,rho,T,A)
% % Initial conditions
% I1=((L/6)/dx)+1;
% I2=((L/2)/dx)+1;
% for i=1:1:I1
%     rho(i)=1;
%     T(i)=1;
% end
% for i=I1+1:1:I2
%     rho(i)=1-0.366*(x(i)-0.5);
%     T(i)=1-0.167*(x(i)-0.5);
% end
% for i=I2+1:1:N
%     rho(i)=0.634-0.3879*(x(i)-1.5);
%     T(i)=0.833-0.3507*(x(i)-1.5);
% end
v=0.59./(rho.*A);
% Initial profiles for U1, U2 and U3
U1=rho.*A;
U2=0.59*ones(1,N);
U3=U1.*(T./(gamma-1)+0.5*gamma*v.^2);
% Initial profiles for F1, F2 and F3
F1=U2;
F2=((U2.^2)./U1)+((gamma-1)/gamma)*(U3-0.5*gamma*((U2.^2)./U1));
F3=gamma*((U2.*U3)./U1)-0.5*gamma*(gamma-1)*((U2.^3)./(U1.^2));
% Initial mass flow rate
ct=1;
dmdt(:,1)=rho.*A.*v;

% Time steps
nt=1400;

% Outer time loop
for k=1:1:nt
    % Calculation of time step using CFL criterion
    for i=2:1:N-1
        time_step(i-1)=C*dx/(T(i)^0.5+v(i));
    end
    dt=min(time_step);
    rho_old=rho;
    v_old=v;
    T_old=T;
    U1_old=U1;
    U2_old=U2;
    U3_old=U3;
    % Predictor step
    for j=2:1:N-1
        %J2=((gamma-1)/gamma)*(U3(j)-0.5*gamma*U2(j)^2/U1(j))*(log(A(j+1))-log(A(j)))/dx;
        J2=(1/gamma)*rho(j)*T(j)*(A(j+1)-A(j))/dx;
        dU1dt_p(j)=-(F1(j+1)-F1(j))/dx;
        dU2dt_p(j)=-(F2(j+1)-F2(j))/dx+J2;
        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;
        F1(j)=U2(j);
        F2(j)=(U2(j)^2/U1(j))+((gamma-1)/gamma)*(U3(j)-0.5*gamma*U2(j)^2/U1(j));
        F3(j)=gamma*U2(j)*U3(j)/U1(j)-0.5*gamma*(gamma-1)*U2(j)^3/U1(j)^2;
    end
    rho=U1./A;
    T=(gamma-1)*((U3./U1)-0.5*gamma*v.^2);
    % Corrector step
    for j=2:1:N-1
        %J2=((gamma-1)/gamma)*(U3(j)-0.5*gamma*U2(j)^2/U1(j))*(log(A(j))-log(A(j-1)))/dx;
        J2=(1/gamma)*rho(j)*T(j)*(A(j)-A(j-1))/dx;
        dU1dt_c(j)=-(F1(j)-F1(j-1))/dx;
        dU2dt_c(j)=-(F2(j)-F2(j-1))/dx+J2;
        dU3dt_c(j)=-(F3(j)-F3(j-1))/dx;
    end
    % Calculation of the average time derivatives
    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 j=2:1:N-1
        U1(j)=U1_old(j)+dt*dU1dt(j);
        U2(j)=U2_old(j)+dt*dU2dt(j);
        U3(j)=U3_old(j)+dt*dU3dt(j);
        F1(j)=U2(j);
        F2(j)=(U2(j)^2/U1(j))+((gamma-1)/gamma)*(U3(j)-0.5*gamma*U2(j)^2/U1(j));
        F3(j)=gamma*U2(j)*U3(j)/U1(j)-0.5*gamma*(gamma-1)*U2(j)^3/U1(j)^2;
    end
    % Boundary conditions
    % Inlet
    U1(1)=A(1);
    U2(1)=2*U2(2)-U2(3);
    U3(1)=U1(1)*((1/(gamma-1))+0.5*gamma*(U2(1)/U1(1))^2);
    F1(1)=U2(1);
    F2(1)=(U2(1)^2/U1(1))+((gamma-1)/gamma)*(U3(1)-0.5*gamma*U2(1)^2/U1(1));
    F3(1)=gamma*U2(1)*U3(1)/U1(1)-0.5*gamma*(gamma-1)*U2(1)^3/U1(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);
    F1(N)=U2(N);
    F2(N)=(U2(N)^2/U1(N))+((gamma-1)/gamma)*(U3(N)-0.5*gamma*U2(N)^2/U1(N));
    F3(N)=gamma*U2(N)*U3(N)/U1(N)-0.5*gamma*(gamma-1)*U2(N)^3/U1(N)^2;
    % Calculation of the primitive variables
    rho=U1./A;
    v=U2./U1;
    T=(gamma-1)*((U3./U1)-0.5*gamma*v.^2);
    p=rho.*T;
    % Calculation of residuals at the throat
    v_residual(k)=max(abs(v-v_old)/dt);
    rho_residual(k)=max(abs(rho-rho_old)/dt);
    T_residual(k)=max(abs(T-T_old)/dt);
    % Calculation of variables atthe throat section
    rho_throat(k)=rho(0.5*(N+1));
    v_throat(k)=v(0.5*(N+1));
    T_throat(k)=T(0.5*(N+1));
    % Calculation of mass flow rate
    if(k==50||k==100||k==150||k==200||k==700)
        ct=ct+1;
        dmdt(:,ct)=rho.*A.*v;
    end
    % Plotting
    % 1.  Residuals vs time steps
    figure(3)
    semilogy(rho_residual,\'-k\');
    hold on;
    semilogy(v_residual,\'--k\');
    semilogy(T_residual,\'-.k\');
    xlabel(\'Number of time steps\');
    ylabel(\'Residuals\');
    legend(\'|d\\rho/dt|_{av}\',\'|dv/dt|_{av}\',\'|dT/dt|_{av}\');
    % 2. Variables vs time steps
    figure(4)
    plot(rho_throat,\'-k\');
    hold on;
    plot(v_throat,\'--k\');
    plot(T_throat,\'-.k\');
    hold off;
    xlabel(\'Number of time steps\');
    ylabel(\'Non-dimensional variables at the throat\');
    legend(\'\\rho/\\rho_0\',\'v/a_0\',\'T/T_0\');
end
% Output
fprintf(\'Conservative form \\n \\n\');
fprintf(\'Non-dimensional density at the throat is=%f \\n\',rho(0.5*(N+1)));
fprintf(\'Exact analytical value of non-dimensional density is=0.634 \\n\');
fprintf(\'Non-dimensional temperature at the throat is=%f \\n\',T(0.5*(N+1)));
fprintf(\'Exact analytical value of non-dimensional temperature is=0.833 \\n\')
fprintf(\'Non-dimensional pressure at the throat is=%f \\n\',p(0.5*(N+1)));
fprintf(\'Exact analytical value of non-dimensional pressure is=0.528 \\n\');
fprintf(\'Mach number at the throat is=%f \\n\',v(0.5*(N+1))/T(0.5*(N+1))^0.5);
fprintf(\'Exact analytical value of Mach number is=1 \\n \\n\');
end

Quasi one-dimensional analysis of subsonic-supersonic isentropic nozzle flow using MacCormack s method

Ship work like this
in your own course.

Talk to a career counsellor

Ship work like thisin your own course.

Talk to a career counsellor

Ship work like this
in your own course.