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

Contents:

• Motivation for doing the project.
• Governing equations: Conservative and Non-conservative forms
• 1D supersonic nozzle flow
  • Numerical procedure
  • Explaining the Macormack method.
• MATLAB codes and functions with initializing the solutions
• Final results that need to be captured
• Results:
  • Compare numerical results with analytical solution
  • Compare conservative and non conservative forms
  • Which method is faster

Motivation

So far we have learnt some ideas on Governing equations, their forms (w/o knowing their actual means), FDM, Consistency of a Numerical scheme, Discretization (FDM), Time integration, Stability and more or less few. Now we put all of these together.

Governing equations

In Conservative form we consider control volume is staionary whereas in Non-Conservative form CV is moving. Both are mathematically same and can easily converted into another easily and excatlty. They have difference only in numerical discretization.

1D Supersonic nozzle flow:

The convergence(subsonic) section connected with reservoir, here mach number is less than one. At throat area, the thermodynamic properties will become critical property as Mach number equal to 1. The sonic flow at the throat means that the local velocity at this point is equal to the local speed of sound and at divergence section(supersonic) the flow expands isentropically to supersonic speed at the nozzle exit. The corresponding notations are clearly written in the figure itself.

The flow is locally subsonic in the convergence section of the nozzle, sonic at the throat and supersonic at the divergence section.

The governing continuity, momentum and energy equation for the quasi 1D steady isentropic flow can be expressed as

Continuity `rho_1V_1A_1 = rho_2V_2A_2`

Momentum `p_1A1 + rho_1V_1^2 + intpdA = p_2A2 + rho_1V_2^2`

Energy `h_1 + V_1^2/2 = h_2 + V_2^2/2`

The real physics of the problem lies on the integral forms of these equations. We will write the equations in their Non-conservative and conservative forms as per Anderson book unit 7.

Non-conservative form

Continuity `(del(rhoA))/(del t) + (del(rhoAV))/(delx) = 0`

Momentum `rho (delu)/(delt) + rhou(delu)/(delx) = -(delp)/(delx)`

Energy `rhoc_v (delT)/(delt) + rhoVc_v (delT)/(delx) = -p(delV)/(delx) - pV(dellnA)/(delx)`

Conservative form

Contiunity `(delrhoA)/(delt) = - (delrhoAV)/(delx)`

Momentum `(delthoAV)/(delt) = - del/(delx)[rhoAV^2 + (1/gamma)pA] + (1/gamma)p(delA)/(delx)`

Energy `del/(delt)[rho(e/(gamma-1) + gamma/2 V^2)A] = -del/(delx)[rho(e/(gamma-1) + gamma/2 V^2)VA + pAV]`

Using Maccormack method, we can solve both of these sets of equations. In this method we use two rounds of iteration namely predictor and corrector steps which solve the premitive variables using FDS and BDS respectively, then we use arthimetic mean of both results.

Following are the sets of Matlab programs with defined initial and boundary conditions as per Anderson book with separate functions for Non-conservative and conservative terms.

Main Program

% Simulation of quasi one dimenstional convergent divergent nozzle
% using macormack method

clear 
clc


%% Discretization
L = 3;   % Total length of Nozzle
n = 31;   % Number of Grid points
x = linspace(0,L,n);   
dx = x(2)-x(1);   % Grid size
throat = (n+1)/2;
% profile of Nozzle and parameters
a =1+2.2*(x-1.5).^2;  % Area parabolic function for x 

% fluid properties
gama = 1.4;
C = 0.5;    % CFL number
nt = 1400;    % number of iterations

tic
[m_Nc,m_Nc_t, M_Nc, M_Nc_t, rho_Nc, rho_Nc_t, P_Nc, P_Nc_t, T_Nc, T_Nc_t,V_Nc, V_Nc_t] = NonConservative(x,n,dx,gama,nt,a,C,throat);
t_non_conservative = toc;
fprintf(\'time consumed for non conservative form = %d\\n\', t_non_conservative);

tic
[m_c,m_c_t, M_c, M_c_t, rho_c, rho_c_t, P_c, P_c_t, T_c, T_c_t,V_c, V_c_t] = Conservative(x,n,dx,gama,nt,a,C,throat);
t_conservative = toc;
fprintf(\'time consumed for conservative form = %d\', t_conservative);

% %% Grid dependence test
% % Calculating primitive values at throat
terms = {\'Nondimensional density\'; \'Nondimentional temperature\'; \'Nondimentional pressure\'; \'Mach number\'};
exact_solution = [0.634; 0.8333; 0.528; 1];
non_conservative_terms = [rho_Nc(throat); T_Nc(throat); P_Nc(throat); M_Nc(throat)];
conservative_terms = [rho_c(throat); T_c(throat); P_c(throat); M_c(throat)];
table(terms, exact_solution, conservative_terms, non_conservative_terms)

figure(3)
plot(x, m_Nc)
hold on
plot(x,m_c)
title(\'comparison of non-dimentional mass flow rates\')
legend (\'Non_conservative form\', \'conservative form\')
xlabel(\'x/L\')
ylabel(\'Nondimensional mass flow rate\')

figure(4)
subplot(4,1,1)
plot(x,M_Nc)
xlabel(\'x-axis\')
ylabel(\'Mach number\')
subplot(4,1,2)
plot(x,rho_Nc)
xlabel(\'x-axis\')
ylabel(\'Density ratio\')
subplot(4,1,3)
plot(x,P_Nc)
xlabel(\'x-axis\')
ylabel(\'Pressure ratio\')
subplot(4,1,4)
plot(x,T_Nc)
xlabel(\'x-axis\')
ylabel(\'Temperature ratio\')
suptitle(\'Solution for Non conservative form at steady state\')

figure(5)
subplot(4,1,1)
plot(x,M_c)
xlabel(\'x-axis\')
ylabel(\'Mach number\')
subplot(4,1,2)
plot(x,rho_c)
xlabel(\'x-axis\')
ylabel(\'Density ratio\')
subplot(4,1,3)
plot(x,P_c)
xlabel(\'x-axis\')
ylabel(\'Pressure ratio\')
subplot(4,1,4)
plot(x,T_c)
xlabel(\'x-axis\')
ylabel(\'Temperature ratio\')
suptitle(\'Solution for Conservative form at steady state\')

Function for Non_Conservative form

% Non Conservative function

function [m_Nc,m_Nc_t, M_Nc, M_Nc_t, rho_Nc, rho_Nc_t, P_Nc, P_Nc_t, T_Nc, T_Nc_t,V_Nc, V_Nc_t] = NonConservative(x,n,dx,gama,nt,a,C, throat);
% Initial conditions along x
rho_Nc = 1-0.3146.*x;  % Density of fluid as a function x
T_Nc = 1-0.2314.*x;    % Temperature of fluid as a function x
V_Nc = (0.1+1.09*x).*T_Nc.^0.5;  % Velocity of fluid as a function x
throat = 16;
% time step


%     dt(i) = C.*(dx./((T(i)).^0.5 +V(i)));
Dt = min(C*(dx./(T_Nc.^0.5 +V_Nc)));
for k = 1:1400
    
    % Values at time t
    rho_old = rho_Nc;
    T_old = T_Nc;
    V_old = V_Nc;
    %a_old = a;
    
    
    
    % Maccormack method
    for j = 2:n-1  % Loop creation for internal grid points
        dvdx = (V_Nc(j+1)-V_Nc(j))/dx;
        dlog_a_dx = (log(a(j+1))-log(a(j)))/dx;
        drhodx = (rho_Nc(j+1)-rho_Nc(j))/dx;
        dTdx = (T_Nc(j+1)-T_Nc(j))/dx;
        
        % Contuinity equation
        drhodt_p(j) = -(rho_Nc(j)*dvdx)-(rho_Nc(j)*V_Nc(j)*dlog_a_dx)-(V_Nc(j)*drhodx);
        % Momentum equation
        dvdt_p(j) = -V_Nc(j)*dvdx-(1/gama)*(dTdx+(T_Nc(j)/rho_Nc(j))*(drhodx));
        % Energy equation
        dTdt_p(j) = -V_Nc(j)*dTdx -(gama-1)*T_Nc(j)*(dvdx+(V_Nc(j)*dlog_a_dx));
        
        % Solution update at time step \'t+dt\'
        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 step
    for j = 2:n-1   % Loop creation for internal grid points
        dvdx = (V_Nc(j)-V_Nc(j-1))/dx;
        dlog_a_dx = (log(a(j))-log(a(j-1)))/dx;
        drhodx = (rho_Nc(j)-rho_Nc(j-1))/dx;
        dTdx = (T_Nc(j)-T_Nc(j-1))/dx;
        
        % Contuinity equation
        drhodt_c(j) = -(rho_Nc(j)*dvdx)-(rho_Nc(j)*V_Nc(j)*dlog_a_dx)-(V_Nc(j)*drhodx);
        % Momentum equation
        dvdt_c(j) = -V_Nc(j)*dvdx-(1/gama)*(dTdx+(T_Nc(j)/rho_Nc(j))*(drhodx));
        % Energy equation
        dTdt_c(j) = -V_Nc(j)*dTdx -(gama-1)*T_Nc(j)*(dvdx+(V_Nc(j)*dlog_a_dx));
    end
    % Arithmetic Mean of pde
    drhodt_avg =0.5*(drhodt_p + drhodt_c);
    dvdt_avg = 0.5*(dvdt_p + dvdt_c);
    dTdt_avg = 0.5*(dTdt_p + dTdt_c);
    % Final solution update
    for i = 2:n-1
        rho_Nc(i) = rho_old(i) +drhodt_avg(i)*Dt;
        V_Nc(i) = V_old(i) +dvdt_avg(i)*Dt;
        T_Nc(i) = T_old(i) +dTdt_avg(i)*Dt;
    end
    
    
    %BC
    V_Nc(1) = 2*V_Nc(2)-V_Nc(3);
    V_Nc(n) = 2*V_Nc(n-1)-V_Nc(n-2);
    rho_Nc(n) = 2*rho_Nc(n-1)-rho_Nc(n-2);
    T_Nc(n) = 2*T_Nc(n-1)-T_Nc(n-2);
    % Calculating
    m_Nc = rho_Nc.*V_Nc.*a;
    M_Nc = V_Nc./(T_Nc.^0.5);
    P_Nc = rho_Nc.*T_Nc;
    
   % flow parameter defining
    P_Nc_t(k) = P_Nc(throat);
    T_Nc_t(k) = T_Nc(throat);
    V_Nc_t(k) = V_Nc(throat);
    rho_Nc_t(k) = rho_Nc(throat);
    m_Nc_t(k) = m_Nc(throat);
    M_Nc_t(k) = M_Nc(throat);
    
    %comparision of Mass Flow rates at different times
    figure(1)
    hold on
    if k ==50
        plot(x,m_Nc,\'y\', \'linewidth\', 1.5)
    elseif k ==100
        plot(x,m_Nc,\'g\', \'linewidth\', 1.5)
    elseif k == 200
        plot(x,m_Nc,\'b\', \'linewidth\', 1.5)
    elseif k ==700
        plot(x,m_Nc,\'--\', \'linewidth\', 1.5)
    end
    title (\'Variation of mass flow distribution in non-conservative form\');
    xlabel(\'non-dimentional length of the nozzle\');
    ylabel(\'mass flow rate\');
    grid on;
end
end

Function for Conservative form

function [m_c,m_c_t, M_c, M_c_t, rho_c, rho_c_t, P_c, P_c_t, T_c, T_c_t,V_c, V_c_t] = Conservative(x,n,dx,gama,nt,a,C, throat);


for j = 1:n
    if (x(j)>=0 && x(j)<= 0.5)
        rho_c(j) =1;
        T_c(j) = 1;
    elseif (x(j)>=0.5 && x(j)<= 1.5)
        rho_c(j) =1-(0.366*(x(j)-0.5));
        T_c(j) = 1-(0.167*(x(j)-0.5));
    elseif (x(j)>=1.5 && x(j)<= 3.5)
        rho_c(j) =0.634-(0.3879*(x(j)-1.5));
        T_c(j) = 0.833-(0.3507*(x(j)-1.5));
    end
end

V_c = 0.59./(rho_c.*a);
m_c = rho_c.*a.*V_c;
u1 = rho_c.*a;
u2 = u1.*V_c;
u3 = u1.*((T_c/(gama-1))+(gama/2).*(V_c.^2));
% dt = 0.002; %min(C.*(dx./((T.^0.5)+V)));  % time step


for k = 1:nt
    Dt = T_c.^0.5;
    dt = min(C*(dx./(Dt+V_c)));
    u1_old = u1;
    u2_old = u2;
    u3_old = u3;
    
    
    % flux vector
    f1 = u2;
    f2_1 = (u2.^2)./u1;
    f2_2 = u3-((gama/2)*f2_1);
    f2 = f2_1 +((gama-1)/gama)*f2_2;
    f3_1 = (u2.*u3)./u1;
    f3_2 = (u2.^3)./(u1.^2);
    f3 = gama.*f3_1 -((gama*(gama-1))/2).*f3_2;
    
    % predictor step
    for i = 2:n-1
        J2_p(i) = (1/gama)*(rho_c(i)*T_c(i))*((a(i+1)-a(i))/dx);
        du1dt_p(i) = -(f1(i+1)-f1(i))/dx;
        du2dt_p(i) = (-(f2(i+1)-f2(i))/dx) + J2_p(i);
        du3dt_p(i) = -(f3(i+1)-f3(i))/dx;
        
        %Solution update
        u1(i) = u1(i) + du1dt_p(i)*dt;
        u2(i) = u2(i) + du2dt_p(i)*dt;
        u3(i) = u3(i) + du3dt_p(i)*dt;
        
    end
    
    % Updating values
    rho_c = u1./a;
    V_c = u2./u1;
    T_c = (gama-1)*((u3./u1)-(gama/2)*(V_c.^2));

    
    f1 = u2;
    f2_1 = (u2.^2)./u1;
    f2_2 = u3-((gama/2)*f2_1);
    f2 = f2_1 +((gama-1)/gama)*f2_2;
    f3_1 = (u2.*u3)./u1;
    f3_2 = (u2.^3)./(u1.^2);
    f3 = gama.*f3_1 -((gama*(gama-1))/2).*f3_2;
    
    
    
    
    % Correcter step
    for i = 2:n-1
        
        J2_c(i) = (1/gama).*(rho_c(i).*T_c(i)).*((a(i)-a(i-1))/dx);
        
        du1dt_c(i) = -(f1(i)-f1(i-1))/dx;
        du2dt_c(i) = (-(f2(i)-f2(i-1))/dx) + J2_c(i);
        du3dt_c(i) = -(f3(i)-f3(i-1))/dx;
    end
    
    % Average values
    du1dt_avg = 0.5*(du1dt_p + du1dt_c);
    du2dt_avg = 0.5*(du2dt_p + du2dt_c);
    du3dt_avg = 0.5*(du3dt_p + du3dt_c);
    
    for z = 2:n-1
        u1(z) = u1_old(z) + du1dt_avg(z)*dt;
        u2(z) = u2_old(z) + du2dt_avg(z)*dt;
        u3(z) = u3_old(z) + du3dt_avg(z)*dt;
    end
    
    % Inlet Boundary conditons
    u1(1) = a(1);
    u2(1) = 2*u2(2) - u2(3);
    u3(1) = u1(1)*((T_c(1)/(gama-1)) + (gama/2)*(V_c(1)^2));
    
    % Outlet Boundary conditions
    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_c = u1./a;
    V_c = u2./u1;
    T_c = (gama-1)*((u3./u1)-(gama/2)*(V_c.^2));
    
    M_c = V_c./(T_c.^0.5);
    m_c = u2;
    P_c = rho_c.*T_c;
   
    
    rho_c_t(k) = rho_c(throat);
    M_c_t(k) = M_c(throat);
    P_c_t(k) = P_c(throat);
    T_c_t(k) = T_c(throat);
    V_c_t(k) = V_c(throat);
    m_c_t(k) = m_c(throat);
    
    figure(2)
    hold on
    if k ==0
        plot(x,m_c,\'--\', \'linewidth\', 1.5)
    elseif k ==100
        plot(x,m_c,\'r\', \'linewidth\', 1.5)
    elseif k == 200
         plot(x,m_c,\'k\', \'linewidth\', 1.5)
    elseif k ==300
        plot(x,m_c,\'c\', \'linewidth\', 1.5)
    elseif k ==700
        plot(x,m_c,\'b\', \'linewidth\', 1.5)

    end
    title (\'Variation of mass flow distribution in conservative form\');
    xlabel(\'non-dimentional length of the nozzle\');
    ylabel(\'mass flow rate\');
    grid on;
end
end

Plots and Observations

In the above two plots, Variation of mass flow distribution through nozzle at different times duting the time marching process for Non-conservative and conservative form of governing equations.

Let us compare the steady-state vatiations of mass flow obtained with the nonconservative and conservative forms of the governing equations (both after 1400 time steps- well beyond the time required to converge to the steady state). This comparision shows that the scale of the ordinate for mass flow is greatly magnified. Hence we see that the steady state mass flow distribution predicted by the conservation form of the equations is much more satisfactory than that obtained from the nonconservation form, which is because

• The conservation form gives a distribution that is much closer to being a constant. In contrast, the nonconservative results have 9on the magnified scale) a sizable variation, with some spurious oscillations at both inflow and outflow boundaries. Of course, on the practical basis, these variations are not apparent and mass flow essentially appears to be constant.
• The steady state mass flow results obtained with the conservative form are, on the whole, much closer to the exact analytical solution of `rho\'A\'V\' = 0.579`

Figure shows a detailed comparison of the steady state mass flow variation(on a magnified scale) obtained with the nonconservative and conservative forms fo the governing flow equations.

In the last two plots, it has been observed that the Mach number increases across the nozzle and pressure, temperature, density ratios are decreases across the non dimentional length of nozzle.

Grid Dependence test:

Note that the noncoservative form results are distinctly closer to the exact values. However, by doubling the number of grid points, the steady state numericals results in conservative form of governing equations are slightly closer to the exact analytical values (but still not as close as the nonconservative form with half as many grid points). For this projects, we have grid indepence with 31 grid points.

Results from the project

For the given flow problem, we cannot establish a clear superiority of the conservation form of the governing equations over the non- conservation from. However, we can make following observations:

1. The conservation form yeilds a better mass flow distribution. The conservation from simply does a better job of conserving mass.
2. The non-consevartive form leads to smallar residuals. The amount by which the computational time is slighly less, however it merely depends on the numerical algorithms in the first place.
3. There is no clear indication that which form yeilds better accuracy, though nonconservative form seems to produce slightly more accurate results for the primitive variables, and the conservative form seems to produce slightly more accurate results for flux variables. This results in either case are certainly satisfacotry.

Which one is faster?

Comparing the amount of calculation effort to achieve a solution, we can note that the solution is conservative form requires marginally more work. Most of this due to the need to decode the primitive variables from the flux variables: such decoding is not necessary when we are solving nonconservative form.