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

Simulating quasi 1-D supersonic flow through a converging diverging nozzle using MacCormack method

Objective:

• Deducing the conservative and non-conservative forms of governing equations and implementing MacCormack method to solve the equations and comparing the results.
• Computing the steady state distribution of fluid variables in the nozzle and the time variation of the variables.
• Studying the variation of mass flow rate inside the nozzle at different time steps under influence of time marching.

Theory:

A nozzle is a type of device designed to control the direction or characteristics of fluid flow at the exit of a pipe (or such system). A nozzle is often a pipe or tube of varying cross-sectional area, thus can be used to direct or modify flow of a fluid.

A convergent-divergent nozzle (also known as CD nozzle or de Laval nozzle) is a device used to accelerate a hot, pressurized gas passing through it to a higher supersonic speed in the axial direction by converting heat energy of the flow to kinetic energy. Hence, such nozzles are widely used in steam turbines, supersnic jet engines and rocket engine nozzles. Its operation relies on the different properties of gases flowing at subsonic, sonic and supersonic speeds.

In the field of CFD, MacCormack method is a prominent discretization scheme for obtaining numerical solutions of hyperbolic partial differential equations. (based on second-order finite difference method). It is the variation of two-step Lax-Wendroff scheme but much simple in application.

In this project, the problem is taken from JD Anderson's book on CFD.

The symmetrical area function for the CD nozzle is given as,

`(A/(A∗))^2=1+2.2(x−1.5)^2 text(@) 0≤x≤3`

where, A∗">A∗A∗ is the area at the throat of the nozzle and x">xx is the distance along the length of the nozzle.

The non-dimensionalized form of the variables required for this problem are,

`x'=x/L; V'=V/a_0; t'=t/(L/a_0); rho'=rho/rho_0; T'=T/T_0; p'=p/p_0; A'=A/(A∗)`

Governing Equations:

1. Non-Conservative form

The non-conservative form of the governing equation which represent continuity, momentum and energy in terms of non-dimensionless variables,

Continuity: `(del rho')/(del t')=−rho' (del V')/(del x')−rho 'V'(del (lnA'))/(del x')−V'(del rho')/(del x')`

Momentum: `(delV')/(del t')=−V'(del V')/(del x')−1/gamma[(delT')/(del x')+(T')/(rho')(del rho')/(del x')]`

Energy: `(del T')/(del t')=−V'(del T')/(del x')−(gamma−1)T'[(del V')/(del x')+V'(del(lnA'))/(del x')]`

The initial conditions as provided in problem statement were,

T′(t=0)=1-0.2314x′">`rho' (t = 0) = 1- 0.3146x'`

T′(t=0)=1-0.2314x′">`T'(t = 0) = 1-0.2314x'`

T′(t=0)=1-0.2314x′">`V'(t= 0) = (0.1 + 1.09x') sqrt(T'=(t = 0))`

2. Conservative Form

The conservative form of governing equation which represents continuity, momentum and energy in non-dimensionalized form is given as,

Continuity: `(del (rho'A'))/(del t')+(del (rho'A'V'))/(del x')=0`

Momentum: `(del (rho'A'V'))/(del t')+(del(rho'A'V'^2)+(1/gamma)rho'A')/(del x')=1/gamma p'(del A')/(del x')`

Energy: `(del [rho′ ((T′)/(gamma-1))+((gamma)/2)V′^2)A′]/(delt′)+(del[rho′((T′)/(gamma-1)+((gamma)/2)V′^2)V′A′+rho′A′V′])/(delx′)=0`

In order to simplify these equations, we substitute the following terms, which include the solution vectors `U_i`, the flux vectors and the source term `J_2`.

`U_1=rho′A′`
`U_2=rho′A′V′`
`U_3=rho′[(T′)/(gamma-1)+gamma/2V′^2]A′`
`F_1=rho′A′V′=U_2`
`F_2=rho′A′V′^2+1/γrho′A′=U_2^2/U_1+(gamma-1)/gamma[U_3-gamma/2U_2^2/U_1]`
`F3=rho′((T′)/(gamma-1)+gamma/2V′^2)V′A′+rho′A′V′=gamma(U_2U_3)/U_1-(gamma(gamma-1))/2U_2^3/U_1^2`
`J2=1/gamma rho′.(delA′) /(del x′)`

After substituting these terms in governing equations, we have,

`(del U_1)/(delt′)=-(delF_1)/(del x′)`
`(del U_2)/(delt′)=-(delF_2)/(del x′)+ J_2`
`(del U_3)/(delt′)=-(delF_3)/(del x′)`

The flow variables are calculated on numerical solution of the above terms as,

`rho′=U_1/(A′)`
`V′=U_2/U_1`
`T′=(gamma-1)[U_3/U_1-gamma/2V′^2]`
`p′=rho′T′`

The initial conditions as provided in the problem statement are,

`rho′(t=0)={(1.0, for ,0, ≤, x′,< ,0.5), (1.0-0.366(x′-0.5), for, 0.5, ≤, x′,< ,1.5),(0.634-0.3879(x′-1.5),for,1.5,≤,x′,≤,3):}`

`T′(t=0)={(1.0 , for ,0, ≤ , x′ ,< ,0.5),(1.0-0.167(x′-0.5),for,0.5,≤,x′,<,1.5),(0.833-0.3507(x′-1.5),for,1.5,≤,x′,≤,3):}`

`V′(t=0)=U_2/(rho′(t=0)A′)=0.59/(rho′(t=0)A′)` (considering a constant mass flow rate of `U_2=rho′A′V′=0.59` just for the initial condition)

MacCormack method of solving:

The MacCormack numerical scheme is a common method to solve flow equations. It basically consists of four steps:

1. Predictor step
2. Intermediate solution update
3. Corrector step
4. Final solution update

In the predictor step, the spatial derivative are calculated using the forward difference scheme. This time derivative from the governing equations obtained is recorded as delu(delt)p">delu(delt)pdelu(delt)p. Consider, the continuity equation for the non-conservative form. From the predictor step, we obtain,

`((del rho_i)/(del t))_P=-rho′_i[(V′_(i+1)-V′_i)/(Deltax′)]-rho′_iV′_i[(lnA′_(i+1)-ln(A′_i))/(Deltax′)]-V′_i[(rho′_(i+1)-rho′_i)/(Deltax′)]`

During the intermediate solution update, the variable is updated based on the time derivative obtained from the predictor step as:

`rho_(i,P)^(t+Delta t)=rho′_i+((del rho′_i)/(del t′))_P Delta t`

Using these updated values, the governing equations are solved again using backward differencing in space in the corrector step. The time-derivative obtained during the corrector step is stored as `((del u)/(del t))_C`.

`((del rho_i)/(delt))_C=-rho′_(i,P)[(V′_(i,P)-V′_(i-1,P))/(Delta x′)]-rho′_(i,P)V′_(i,P)[(lnA′_(i,P)-ln(A′_(i-1,P)))/(Deltax′)]-V′_(i,P)[(rho′_(i,P)-rho′_(i-1,P))/(Deltax′)]`

For the final solution update, we use the average of the time-derivatives obtained using the predictor and the corrector steps. Since we are using forward difference in predictor step and backward difference in corrector step, the average of the time-derivative for the final update, the overall scheme becomes a second order numerical scheme.

`((del rho)/(delt))_(av)=1/2[((del rho)/(del t))_P+((delrho)/(delt))_C]`

`rho_i^(t+Delta t)=rho_i^t+((del rho)/(del t))_(av)Delta t`

Matlab Code:

close all
clear all
clc

% Initializing the domain
L = 3;
nx = 31;
x = linspace(0,L,nx)';
dx = x(2)-x(1);
C = 0.5;
nt = 1400;
gamma = 1.4;

% Nozzle area profile
A = 1 + 2.2*(x-1.5).^2;

% Computing the flow variable through conservative and non-conservative form
[rho1, V1, T1, p1, M1, m1, rho_t1, V_t1, p_t1, T_t1, M_t1, residual1] = nonconservative(x, dx, A, C, nt, gamma);
[rho2, V2, T2, p2, M2, m2, rho_t2, V_t2, p_t2, T_t2, M_t2, residual2] = conservative(x, dx, A, C, nt, gamma);

% Plotting the results
figure(1)
subplot(4,1,1)
plot(x,M1)
grid on; ylim([0 4.5]); yticks(0:1.5:4.5);
ylabel('Mach Number')
title('Flow parameter (non-conservative form of governing equation)')
subplot(4,1,2)
plot(x,p1)
grid on; yticks(0:0.2:1);
ylabel('Pressure Ratio')
subplot(4,1,3)
plot(x,rho1)
grid on; yticks(0:0.2:1);
ylabel('Density Ratio')
subplot(4,1,4)
plot(x,T1)
grid on; ylim([0 1]); yticks(0:0.2:1);
ylabel('Temperature Ratio')
xlabel('X/L')

figure(2)
subplot(4,1,1)
plot(x,M2)
grid on; ylim([0 4.5]); yticks(0:1.5:4.5);
ylabel('Mach Number')
title('Flow parameters (conservative form of governing equation)')
subplot(4,1,2)
plot(x,p2)
grid on; yticks(0:0.2:1);
ylabel('Pressure Ratio')
subplot(4,1,3)
plot(x,rho2)
grid on; yticks(0:0.2:1);
ylabel('Density Ratio')
subplot(4,1,4)
plot(x,T2)
grid on; ylim([0 1]); yticks(0:0.2:1);
ylabel('Temperature Ratio')
xlabel('X/L')

figure(3)
subplot(4,1,1)
plot(1:nt,rho_t1)
hold on
plot(1:nt,rho_t2)
grid on; ylim([0 1]); yticks(0:0.2:1);
ylabel('rho/rho_0')
title('Flow parameter variation with time at the nozzle throat');
legend('Non-conservative','Conservative');
subplot(4,1,2)
plot(1:nt,T_t1)
hold on
plot(1:nt,T_t2)
grid on; yticks(0:0.2:1);
ylabel('T/T_0')
legend('Non-conservative','Conservative');
subplot(4,1,3)
plot(1:nt,p_t1)
hold on
plot(1:nt,p_t2)
grid on; yticks(0:0.2:1);
ylabel('p/p_0')
legend('Non-conservative','Conservative')
subplot(4,1,4)
plot(1:nt,M_t1)
hold on
plot(1:nt,M_t2)
grid on; yticks(0:0.2:2);
ylabel('Mach No.')
xlabel('Time Steps')
legend('Non-conservative','Conservative')

figure(4)
subplot(2,1,1)
plot(1:nt,residual1)
grid on
xlabel('Time Steps')
ylabel('Residuals')
legend('|drho/dt|_{av}','|dV/dt|_{av}','|dT/dt|_{av}')
title('Residuals: Non-conservative form')
subplot(2,1,2)
plot(1:nt,residual2)
grid on
xlabel('Time Steps')
ylabel('Residuals'); yticks(0:0.5:4);
legend('|dU_1/dt|_{av}','|dU_2/dt|_{av}','|dU_3/dt|_{av}')
title('Residuals: Conservative form')

figure(5)
hold on; grid on;
plot(x,m1(50,:),'-x')
plot(x,m1(100,:),'-d')
plot(x,m1(150,:),'-s')
plot(x,m1(200,:),'-*')
plot(x,m1(700,:),'-o')
xlabel('X/L')
ylabel('rhoAV/rho_0A_0V_0')
legend('t = 50/\t','t = 100/\t','t = 150/\t','t = 200/\t','t = 700/\t')
title('Variation of non-dimensional mass flow rate (Non-conservative form)')

figure(6)
hold on; grid on;
plot(x,m2(100,:),'-x')
plot(x,m2(150,:),'-d')
plot(x,m2(200,:),'-*')
plot(x,m2(700,:),'-o')
xlabel('X/L')
ylabel('rhoAV/rho_0A_0V_0')
legend('t = 100/\t','t = 150/\t','t = 200/\t','t = 700/\t')
title('Variation of non-dimensional mass flow rate (Conservative form)')

figure(7)
hold on; grid on;
plot(x,m1(end,:),'LineWidth',3)
plot(x,m2(end,:),'LineWidth',3)
plot(x,0.579*ones(length(x),1),'--k')
xlabel('X/L')
ylabel('rhoAV/rho_0A_0V_0')
legend('Non-conservative','Conservative','Exact')
title('Variation of non-dimensional mass flow rate (Steady State)')

%% Function Statement 

% Solver for Conservative Form of Governing Equation
function [rho, V, T, p, M, m, rho_t, V_t, p_t, T_t, M_t, residual] = conservative(x, dx, A, C, nt, gamma)
    % Initial conditions
    rho_0 = zeros(length(x),1);
    T_0 = rho_0;
    V_0 = T_0;
    for i=1:length(x)
        if x(i)>=0 && x(i)<0.5
            rho_0(i) = 1;
            T_0(i) = 1;
        elseif x(i)>=0.5 && x(i)<1.5
            rho_0(i) = 1 - 0.366*(x(i)-0.5);
            T_0(i) = 1 - 0.167*(x(i)-0.5);
        elseif x(i)>=1.5 && x(i)<=3
            rho_0(i) = 0.634 - 0.3879*(x(i)-1.5);
            T_0(i) = 0.833 - 0.3507*(x(i)-1.5);
        end
        V_0(i) = 0.59/(rho_0(i)*A(i));
    end
    throat = find(A==1);
    rho = rho_0;
    T = T_0;
    V = V_0;
    m = zeros(nt,length(x));
    U1 = rho.*A;
    U2 = rho.*A.*V;
    U3 = rho.*(T./(gamma-1) + 0.5*gamma*V.^2).*A;
    rho_t = zeros(nt,1);
    T_t = rho_t;
    V_t = T_t;
    p_t = V_t;
    M_t = p_t;
    residual = zeros(nt,3);
    for k=1:nt
        U1_old = U1; U2_old = U2; U3_old = U3;
        dt = min(C*dx./(T.^0.5 + V));
        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);
        % Predictor loop
        for i=2:length(x)-1
            dF1dx = (F1(i+1) - F1(i))/dx;
            dF2dx = (F2(i+1) - F2(i))/dx;
            dF3dx = (F3(i+1) - F3(i))/dx;
            dAdx = (A(i+1) - A(i))/dx;
            J2 = (1/gamma)*rho(i)*T(i)*dAdx;

            dU1dt_p(i) = -dF1dx;
            dU2dt_p(i) = -dF2dx + J2;
            dU3dt_p(i) = -dF3dx;

            % Intermediate 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
        rho = U1./A;
        V = U2./U1;
        T = (gamma-1)*((U3./U1) - 0.5*gamma*V.^2);
        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);
        % Corrector loop
        for i=2:length(x)-1
            dF1dx = (F1(i) - F1(i-1))/dx;
            dF2dx = (F2(i) - F2(i-1))/dx;
            dF3dx = (F3(i) - F3(i-1))/dx;
            dAdx = (A(i) - A(i-1))/dx;
            J2 = (1/gamma)*rho(i)*T(i)*dAdx;

            dU1dt_c(i) = -dF1dx;
            dU2dt_c(i) = -dF2dx + J2;
            dU3dt_c(i) = -dF3dx;
        end

        % Average of predictor and corrector step
        dU1dt = 0.5*(dU1dt_p + dU1dt_c);
        dU2dt = 0.5*(dU2dt_p + dU2dt_c);
        dU3dt = 0.5*(dU3dt_p + dU3dt_c);
        residual(k,1) = max(abs(dU1dt));
        residual(k,2) = max(abs(dU2dt));
        residual(k,3) = max(abs(dU3dt));

        % Final solution update
        for i=2:length(x)-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 conditions
        U1(1) = rho(1)*A(1);
        U1(end) = 2*U1(end-1) - U1(end-2);
        U2(1) = 2*U2(2) - U2(3);
        U2(end) = 2*U2(end-1) - U2(end-2);
        U3(1) = U1(1)*((T(1)/(gamma-1)) + 0.5*gamma*V(1)^2);
        U3(end) = 2*U3(end-1) - U3(end-2);
        rho = U1./A;
        V = U2./U1;
        T = (gamma-1)*((U3./U1) - 0.5*gamma*V.^2);
        M = V./(T.^0.5);
        p = rho.*T;
        m(k,:) = rho.*A.*V;
        rho_t(k) = rho(throat); V_t(k) = V(throat); T_t(k) = T(throat);
        p_t(k) = p(throat); M_t(k) = M(throat);
    end
end

%% Function Statement 
% Solver for Non-Conservative Form of Governing Equation
function [rho, V, T, p, M, m, rho_t, V_t, p_t, T_t, M_t, residual] = nonconservative(x, dx, A, C, nt, gamma)
    % Initial conditions
    rho_0 = 1 - 0.3146*x;
    T_0 = 1 - 0.2314*x;
    V_0 = (0.1 + 1.09*x).*T_0.^0.5;
    throat = find(A==1);
    rho = rho_0; T = T_0; V = V_0; m = zeros(nt,length(x));
    rho_t = zeros(nt,1); T_t = rho_t; V_t = T_t; p_t = V_t; M_t = p_t;
    residual = zeros(nt,3);
    for k=1:nt
        rho_old = rho; T_old = T; V_old = V;
        dt = min(C*dx./(T.^0.5 + V));
        % Predictor loop
        for i=2:length(x)-1
            dVdx = (V(i+1) - V(i))/dx;
            dlogAdx = (log(A(i+1))-log(A(i)))/dx;
            drhodx = (rho(i+1) - rho(i))/dx;
            dTdx = (T(i+1) - T(i))/dx;

            drhodt_p(i) = -rho(i)*dVdx - rho(i)*V(i)*dlogAdx - V(i)*drhodx;
            dVdt_p(i) = -V(i)*dVdx - (1/gamma)*(dTdx + (T(i)/rho(i))*drhodx);
            dTdt_p(i) = -V(i)*dTdx - (gamma-1)*T(i)*(dVdx + V(i)*dlogAdx);

            % Intermediate update
            rho(i) = rho(i) + drhodt_p(i)*dt;
            V(i) = V(i) + dVdt_p(i)*dt;
            T(i) = T(i) + dTdt_p(i)*dt;
        end
        % Corrector loop
        for i=2:length(x)-1
            dVdx = (V(i) - V(i-1))/dx;
            dlogAdx = (log(A(i))-log(A(i-1)))/dx;
            drhodx = (rho(i) - rho(i-1))/dx;
            dTdx = (T(i) - T(i-1))/dx;

            drhodt_c(i) = -rho(i)*dVdx - rho(i)*V(i)*dlogAdx - V(i)*drhodx;
            dVdt_c(i) = -V(i)*dVdx - (1/gamma)*(dTdx + (T(i)/rho(i))*drhodx);
            dTdt_c(i) = -V(i)*dTdx - (gamma-1)*T(i)*(dVdx + V(i)*dlogAdx);
        end

        % Average of predictor and corrector step
        drhodt = 0.5*(drhodt_p + drhodt_c);
        dVdt = 0.5*(dVdt_p + dVdt_c);
        dTdt = 0.5*(dTdt_p + dTdt_c);
        residual(k,1) = max(abs(drhodt));
        residual(k,2) = max(abs(dVdt));
        residual(k,3) = max(abs(dTdt));

        % Final solution update
        for i=2:length(x)-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

        % Boundary conditions
        rho(end) = 2*rho(end-1)-rho(end-2);
        V(1) = 2*V(2)-V(3); V(end) = 2*V(end-1)-V(end-2);
        T(end) = 2*T(end-1)-T(end-2);
        M = V./(T.^0.5);
        p = rho.*T;
        m(k,:) = rho.*A.*V;
        rho_t(k) = rho(throat); V_t(k) = V(throat); T_t(k) = T(throat);
        p_t(k) = p(throat); M_t(k) = M(throat);
    end
end

Plots:

Mass Flow nx = 31 ;

Mass Flow nx = 71 ;

Result/Conclusion:

• For comparing the conservative and non-conservative form, the values of the flow variables at the throat of the nozzle are plotted at each time step. It can be observed that at smller time-steps, the behaviour of non-conservative form shows unstable results. (it shows oscillations) but with increase in time step, the difference in conservative and non-conservative approach diminshes.
• If we notice the residuals (average time derivatives of variables) of both the approaches the conservative form converges to a steady state faster than non-conservative form.
• There is an appreciable difference in how the mass flow rate varies with time step, but it can be observed that at the time step of `700 Delta t`, the mass flow becomes constant. Moreover, the conservative form produces an almost constant mass flow rate but the non-conservative approach gave a highly fluctuating mass flow rate.
• On varying the grid size, it is noted that the conservative form almost maintains its consistency but the non conservative form tends to become stable (on increasing grid size).
• Thus MaCormack method provides quite stable and close to exact solutions, hence, it is widely implemented to solve various flow problems.