1D Super-sonic nozzle flow simulation

using Macormack Method

Objective :-

Your objective in this project is to write code solve the 1D supersonic nozzle flow equations using the Macormack Method.

• You will be implementing both the conservative and non-conservative forms of the governing equations
• Write separate functions for conservative and non-conservative forms
• Figure out the minimum number of cycles for which the simulation should be run in order for convergence
• Compare the normalized mass flow rate between the conservative forms and non-conservative forms - make suitable conclusions

MacCormack method

In computational fluid dynamics, the MacCormack method is a widely used discretization scheme for the numerical solution of hyderbolic partial differential equation. This second-order finie difference method was introduced by Robert W. MacCormack in 1969. The MacCormack method is elegant and easy to understand and program.

The MacCormack method is a variation of the two step Lax-Wendroff scheme but is much simpler in application. To illustrate the algorithm, consider the following first order hyperbolic equation

The application of MacCormack method to the above equation proceeds in two steps; a predictor step which is followed by a corrector step.

Predictor step: In the predictor step, a "provisional" value of at time level is estimated as follows

The above equation is obtained by replacing the spatial and temporal derivatives in the previous first order hyperbolic equation using forward differences.

Corrector step: In the corrector step, the predicted value is corrected according to the equation

Note that the corrector step uses backward finite difference approximations for spatial derivative. The time-step used in the corrector step is in contrast to the used in the predictor step.

Replacing the term by the temporal average

to obtain the corrector step as

The MacCormack method is well suited for non-linear equation (Inviscid burgers equation,Euler equation etc.) The order of differencing can be reversed for the time step (i.e., forward/backward followed by backward/forward). For nonlinear equations, this procedure provides the best results. For linear equations, the MacCormack scheme is equivalent to the Lax-Wendroff method.

Unlike first-order upwind scheme, the MacCormack does not introduce diffusive errors in the solution. However, it is known to introduce dispersive errors (Gibbs phenomenon) in the region where the gradient is high.

Main Programe :-

clear all
close all
clc

%Inputs
n = 31;
x = linspace(0,3,n);
dx = x(2)-x(1);
gamma = 1.4;

% Time steps
nt = 1000;
c = 0.5;
tol = 1e-6;

%calling function programe 
[rho1,v1,t1,k1,netmf1,error_rho1,error_v1,error_t1] = non_conserve(n,x,dx,gamma,nt,c,tol);
[rho2,v2,t2,k2,netmf2,error_rho2,error_v2,error_t2] = conserve(n,x,dx,gamma,nt,c,tol);

%Velocity plot
figure(1)
plot(x,v1,'b',x,v2,'g','linewidth',2);
xlabel('Length');
ylabel('Velocity');
legend('Non conservative','Conservative');
title('Non conservative vs Conservative');
grid on 

%Temperature plot
figure(2);
plot(x,t1,'b',x,t2,'g','linewidth',2);
xlabel('Length');
ylabel('Temprature');
legend('Non conservative','Conservative');
title('Non conservative vs Conservative');
grid on

% Density plot
figure(3);
plot(x,rho1,'b',x,rho2,'g','linewidth',2)
xlabel('Length');
ylabel('Density');
legend('Non Conservative','Conservative');
title('Non conservative vs Conservative');
grid on

% Plot for Net mass flow rate 
figure(4)
subplot(2,1,1)
plot(netmf1,'b')
xlabel('Time Step');
ylabel('Net Mass Flow Rate');
title('Non conservative Form');
grid on
subplot(2,1,2)
plot(netmf2, 'g')
xlabel('Time Step');
ylabel('Net Mass Flow Rate');
title('Conservative Form');
grid on


% Error in Velocity with time
figure(5);
subplot(2,1,1)
plot(error_v1,'b')
xlabel('Time Step');
ylabel('Error in velocity');
title('Non Conservative Form');
grid on
subplot(2,1,2)
plot(error_v2,'b')
xlabel('Time Step');
ylabel('Error in velocity');
title('Conservative Form');
grid on


% Error in density with time
figure(6);
subplot(2,1,1)
plot(error_rho1,'b')
axis([0 inf -0.01 0.02])
xlabel('Time Step');
ylabel('Error in Density');
title('Non Conservative Form');
grid on
subplot(2,1,2)
plot(error_rho2,'b')
axis([0 inf -0.005 0.015]);
xlabel('Time step');
ylabel('Error in Density')
title('Conservative Form')
grid on

%Plot for Error in temperature with time
figure(7);
subplot(2,1,1);
plot(error_t1,'b')
axis([0 inf -0.005 0.015]);
xlabel('Time step');
ylabel('Error in Temperature')
title('Non Conservative Form')
grid on
subplot(2,1,2)
plot(error_t2,'b')
axis([0 inf -0.01 0.02]);
xlabel('Time Step');
ylabel('Error in Temperature');
title('Conservative Form');
grid on


    
    

Function programe for Non-conservative :-

%Function for Non-conservative form
function [rho,v,t,k,netmf,error_rho,error_v,error_t] = non_conserve(n,x,dx,gamma,nt,c,tol)

% calculating initial profiles
rho = 1-0.3146*x;
% Temperature
t = 1-0.2314*x;
% Velocity
v = (0.1+0.9*x).*t.^0.5;
% Area of nozzle
a = 1+2.2*(x-1.5).^2; 

% Outer time loop
for k = 1:nt
    
    dt = min(abs(c*dx./(realsqrt(t)+v)));
    
    % saving the old variables
    rho_old = rho;
    t_old = t;
    v_old = v;
    
    % 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;
      
        % continuity method 
        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 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;
      
        % continuity method 
        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
    
    % computate 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 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
    
    % Apply boundary condition
    % Inlet
    v(1) = 2*v(2) - v(3);
    
    %outlet
    v(n) = 2*v(n-1) - v(n-2);
    rho(n) = 2*rho(n-1) - rho(n-2);
    t(n) = 2*t(n-1) - t(n-2);
    
    % Mass flow rate
    inflow = rho(1)*v(1)*a(1);
    outflow = rho(n)*v(n)*a(n);
    netmf(k) = (outflow-inflow);
    
    % Error ion tolerence
    error_rho(k) = max(rho-rho_old);
    error_v(k) = max(v-v_old);
    error_t(k) = max(t-t_old);
    
   
end

        
    
    

Function programe for conservative :-

function[rho,v,t,k,netmf,error_rho,error_v,error_t] = conserve(n,x,dx,gamma,nt,c,tol)

% calculating initial profiles
rho = 1-0.3146*x;
% Temperature
t = 1-0.2314*x;
% Velocity
v = (0.1+0.9*x).*t.^0.5;
% Area of nozzle
a = 1+2.2*(x-1.5).^2; 

% Derivatives
U1 = rho.*a;
U2 = rho.*a.*v;
U3 = rho.*(t/(gamma-1)+(gamma/2)*v.^2).*a;

%outer time loop
for k = 1:nt
    
    % saving old values 
    U1_old = U1;
    U2_old = U2;
    U3_old = U3;
    
    
    % Copying old variable
    rho_old = U1_old./a;
    V_old = U2_old./(rho.*a);
    t_old = (gamma-1)*((U3_old./U1_old)-(gamma/2)*v.^2);
    
    r1 = U2;
    r2 = (U2.^2)./U1+((gamma-1)/gamma)*(U3-(gamma/2)*(U2.^2)./U1);
    r3 = (gamma*U2.*U3)./U1-((gamma*(gamma-1)/2)*(U2.^3)./U1.^2);
    
    for j = 2:n-1
        J2(j) = (1/gamma)*rho(j)*t(j)*(a(j+1)-a(j))/dx;
    end
    
    for i = 1:n
        deltat(i) = (c*(dx/((t(i)^0.5+v(i)))));
    end
    
    dt = min(deltat);
    
    % Predictor Method 
    for j = 2:n-1
        
       % Continuity Equation
       dU1dt_p(j) = -(r1(j+1)-r1(j))/dx;
       
       % Momentum Equation
       dU2dt_p(j) = -(r2(j+1)-r2(j))/dx+J2(j);
       
       % Energy Equation
       dU3dt_p(j) = -(r3(j+1)-r3(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
    
    % Calculating the primary Value
    rho = U1./a;
    v = U2./(rho.*a);
    t = (gamma-1)*((U3./U1)-(gamma/2)*v.^2);
    
     r1 = U2;
    r2 = (U2.^2)./U1+((gamma-1)/gamma)*(U3-(gamma/2)*(U2.^2)./U1);
    r3 = (gamma*U2.*U3)./U1-((gamma*(gamma-1)/2)*(U2.^3)./U1.^2);
    
    for j = 2:n-1
        J2(j) = (1/gamma)*rho(j)*t(j)*(a(j+1)-a(j))/dx;
    end
    
    % Corrector Method 
    for j = 2:n-1
        
       % Continuity Equation
       dU1dt_c(j) = -(r1(j)-r1(j-1))/dx;
       
       % Momentum Equation
       dU2dt_c(j) = -(r2(j)-r2(j-1))/dx+J2(j);
       
       % Energy Equation
       dU3dt_c(j) = -(r3(j)-r3(j-1))/dx;
       
    end
    
    % Computating average time derivative 
    dU1dt = 0.5*(dU1dt_p+dU1dt_c);
    dU2dt = 0.5*(dU2dt_p+dU2dt_c);
    dU3dt = 0.5*(dU3dt_p+dU3dt_c);
    
     % Solution Update 
     for j = 2:n-1
       U1(j) = U1_old(j)+dU1dt(j)*dt;
       U2(j) = U2_old(j)+dU2dt(j)*dt;
       U3(j) = U3_old(j)+dU3dt(j)*dt;
     end
     
     % Applying Boundary Condition
     % Inlet
     U2(1) = 2*U2(2)-U2(3);
          
     % 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);

    % Recalculating Primitive values
    rho = U1./a;
    v = U2./(rho.*a);
    t = (gamma-1)*((U3./U1)-(gamma/2)*v.^2);
    
    % Mass Flow Rate 
    inflow = rho(1)*v(1)*a(1);
    outflow = rho(n)*v(n)*a(n);
    netmf(k) = (outflow-inflow);
    
    % Error tolerence
    error_v(k) = max(v-V_old);
    error_rho(k) = max(rho-rho_old);
    error_t(k) = max(t-t_old);








end

Explanation:-

• The density,temperature and velocity profiles obtained from the method of conservative and non-conservative are good.
• The mass flow distribution of the conservative m,ethod look more stable than the conservative method and its is clear then the mass flow rate of the two different methods are very much different from each other.
• We using non consrvative and conservative formula are se used function loop .
• The function loop iteratiion has been started we plotted the graph so it will print the values.

Output :-

Workspace :-

Conclusion:-

• After comparing both the methods it has been seen that there are a lot of fluctuation in non-conservative flow when comapred to the conservative form.
• The both methods has been seen that there are more fluctuation in non conseravtive flow whem compared to the conservative one.
• The Mass flow distribution is good in conservation form because the mass flow can be seen as dependent variables in the conservative from equations whereas in the case of non-conservative method the mass flow can be considered as a secondary result. hrence mass is more conserved in conservative form.
• The conservatyive form is better in mass flow rate over distribution the length of the nozzle and thus give better mass conservation.
• Non-conservative method takes more time for solving and more number of iterations for convergence.