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

Theoretical Aspects:

The conservative form and non-conservative form are the two type numerical discretization in the computational fluid dynamics. Conservative form comes when we assume the control volume is stationary and on the other hand non-conservative form comes when we assume control volume is moving. Our primary objective for this challenge how can we relate this two method for the simplication and solution generation for the quas 1D subsonic and supersonic nozzle section.

The conservative form and non-conservative form are written in both integral form and PDE form. For the coding the most preferable is partial differential equation form. We are solving the both conservative and non-conservative form using MacCormak's method, which provide the second order accuracy without involvment of second order derivative terms in the equation. Which greatly reduce the mathematical work and coding work for the problem.

In the MacCormak's Method, first we derive the eqaution for the predector step in which we use forward differencing approach and after that we get the solution updates value initial value plus the predector step value multiplied by the time step size. After calculating the solution update term we move the corrector step in which for the same equation we take backward differencing and then take the average of the value, and in the last step we calculate the quantity at next time step by initial size added by the average term multiplied again with step size.

The Governing equation for Non-conservative form is given as:

1) Continuity Equation:

`((delrho')/(delt')) = rho'*((delV')/(delx')) - rho'*V'*((del(LnA))/(delx'))- V'*((delrho')/(delx'))`

2) Momentum Equation:

`((delV')/(delt')) = V'*((delV')/(delx'))-(1/gamma)*(((delT')/(delx'))+((T')/(rho'))*((delrho')/(delx')))`

3) Energy Equation:

`((delT')/(delt')) = V'*((delT')/(delx'))-(gamma-1)*T'*[((delV')/(delx'))+V'*((del(LnA'))/(delx'))]`

The Governing equation for conservative form is given as:

1) Continuity Equation:

`(del(rho'*A'))/(delt')+(del(rho'*A'*V'))/(delx') = 0`

2) Momentum Equation:

`(del(rho'*A'*V'))/(delt') + (del(rho'*A'*V'^2+(1/gamma)*p'*A'))/(delx') = (1/gamma)*p'(delA')/(delx')`

3) Energy Equation:

`(del(rho'*((T')/(gamma-1)+gamma/2*V'^2)A'))/(delt')+(del(rho'*((T')/(gamma-1)+gamma/2*V'^2)*V'*A'+p'*A'*V'))/(delx')=0`

The input for the both form of equation are given and the boundary condition is also provided.

Main Code for Non-Conservative & Conservative Function

clear all
close all
clc

n = 31;                     % No. of gris points
L = 3;                      % Length of Domain
x = linspace(0,L,n);        % Discretization
c = 0.5;                    % Value of courant number
throat = ((n-1)/2)+1;       % Throat dimension
nt = 1400;                  % Number of time steps

% Non-conservative form of Equation
tic
[v,rho,t,m,mach,pr,v_t,rho_t,t_t,m_t,mach_t,pr_t] = non_conservative_equation(n,c,L,nt,throat);
Non_Conservative_time = toc;

tic
[v2,rho2,t2,m2,mach2,pr2,v_t2,rho_t2,t_t2,m_t2,mach_t2,pr_t2] = conservative_equation(n,L,c,nt,throat);
Conservative_time = toc;

%Flow parameters at throat
figure(3)
e1 = subplot(4,1,1);
plot(linspace(1,nt,nt),rho_t);
hold on
plot(linspace(1,nt,nt),rho_t2,'color','r');
hold on
plot(linspace(nt-400,nt+200,600),linspace(0.634,0.634,600),'color','b')
grid on
ylabel("rho'",'Fontsize',12);
h = legend(e1,{'Non-Conservatice form','Conservative form','Analytical solution'});
set(h,'Location','northeastoutside')
title('Timewise Variation of Flow Paramenters At Throat','Fontsize',12,'color','r');


e2 = subplot(4,1,2);
plot(linspace(1,nt,nt),t_t);
hold on
plot(linspace(1,nt,nt),t_t2,'color','r');
hold on
plot(linspace(nt-400,nt+200,600),linspace(0.833,0.833,600),'color','b')
grid on
ylabel("T'",'Fontsize',12);
h = legend(e2,{'Non-Conservatice form','Conservative form','Analytical solution'});
set(h,'Location','northeastoutside')


e3 = subplot(4,1,3);
plot(linspace(1,nt,nt),pr_t);
hold on
plot(linspace(1,nt,nt),pr_t2,'color','r');
hold on
plot(linspace(nt-400,nt+200,600),linspace(0.528,0.528,600),'color','b')
grid on
ylabel("Pr'",'Fontsize',12);
h = legend(e3,{'Non-Conservatice form','Conservative form','Analytical solution'});
set(h,'Location','northeastoutside')


e4 = subplot(4,1,4);
plot(linspace(1,nt,nt),mach_t);
hold on
plot(linspace(1,nt,nt),mach_t2,'color','r');
hold on
plot(linspace(nt-400,nt+200,600),linspace(1.0,1.0,600),'color','b')
grid on
ylabel("Ma'",'Fontsize',12);
h = legend(e4,{'Non-Conservatice form','Conservative form','Analytical solution'});
set(h,'Location','northeastoutside')

%Variation of parameters along x axis
figure(4)
q1 = subplot(4,1,1);
plot(x,mach);
hold on
plot(x,mach2,'*');
grid on;
ylabel("Ma'",'Fontsize',12)
h = legend(q1,{'Non-conservative form','Conservative form'});
set(h,'Location','northeastoutside')
title('Variation of Flow parameters along x axis','Fontsize',12,'color','r')

q2 = subplot(4,1,2);
plot(x,pr);
hold on
plot(x,pr2,'*');
grid on;
ylabel("Pr'",'Fontsize',12)
h = legend(q2,{'Non-conservative form','Conservative form'});
set(h,'Location','northeastoutside')
title('Variation of Flow parameters along x axis','Fontsize',12,'color','r')

q3 = subplot(4,1,3);
plot(x,rho);
hold on
plot(x,rho2,'*');
grid on;
ylabel("rho'",'Fontsize',12)
h = legend(q3,{'Non-conservative form','Conservative form'});
set(h,'Location','northeastoutside')
title('Variation of Flow parameters along x axis','Fontsize',12,'color','r')

q4 = subplot(4,1,4);
plot(x,t);
hold on
plot(x,t2,'*');
grid on;
ylabel("T'",'Fontsize',12)
h = legend(q4,{'Non-conservative form','Conservative form'});
set(h,'Location','northeastoutside')
title('Variation of Flow parameters along x axis','Fontsize',12,'color','r')

% Mass flow Comparison
figure(5)
plot(x,m,'r','linewidth',1.5)
hold on
plot(x,m2,'b','linewidth',1.5)
hold on
plot(x,linspace(0.579,0.579,n),'color','m')
xlabel('Non-Dimensional Distance of Nozzle x','fontsize',12);
ylabel('Non-Dimensional Mass flow rate','fontsize',12);
legend({'Non-conservative form','Conservative form','Analytical Solution'});
title('Non-Dimensional Mass Flow Rate Comparison','color','r','fontsize',12)
axis([0 3 0.575 0.600]);
grid on

The main function code for the Non-conservative form is given as:

function [v,rho,t,m,mach,pr,v_t,rho_t,t_t,m_t,mach_t,pr_t] = non_conservative_equation(n,c,L,nt,throat)

%inputs
x = linspace(0,L,n);
dx = x(2)-x(1);

%Initial conditions
rho = 1-0.3146*x;
t = 1-0.2314*x;
v = (0.1+1.09*x).*t.^0.5;
a = 1+2.2*(x-1.5).^2;
gamma = 1.4;

%Time step value
for i = 1:n
    b(i) = c*(dx/(t(i)^0.5+v(i)));
end
dt = min(abs(b));

% Time Loop
for k = 1:nt
    
    rho_old = rho;
    v_old = v;
    t_old = t;
    
    %Predictor methods
    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; %keep track of this bracket
        
        %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 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; %keep track of this bracket
        
        %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
    %Average derivative
    drhodt = 0.5*(drhodt_p + drhodt_c);
    dvdt = 0.5*(dvdt_p + dvdt_c);
    dtdt = 0.5*(dtdt_p + dtdt_c);
    
    % Final time step
    for p = 2:n-1
        rho(p)= rho_old(p) + drhodt(p)*dt;
        v(p) = v_old(p) + dvdt(p)*dt;
        t(p) = t_old(p) + dtdt(p)*dt;
    end
    
    % Apply boundary condtions
    % Inlet BC
    v(1) = 2*v(2)-v(3);
    rho(1) = 1; % Check this
    t(1) = 1; % Check this also
    
    %Outlet BC
    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);
    
    % Calculation of solutions
    
    m = rho.*v.*a;
    mach = v./((t).^0.5);
    pr = rho.*t;
    
    %throat value
    
    rho_t(k) = rho(throat);
    t_t(k) = t(throat);
    pr_t(k) = pr(throat);
    mach_t(k) = mach(throat);
    v_t(k) = v(throat);
    m_t(k) = m(throat);
    
    %Plotting mass flow rate
    
    figure(1);
    if k == 50
        plot(x, m,'b','linewidth',1.5)
        hold on
    elseif k==100
        plot(x, m,'r','linewidth',1.5)
        hold on
    elseif k==150
        plot(x, m,'k','linewidth',1.5)
        hold on
    elseif k==200
        plot(x, m,'y','linewidth',1.5)
        hold on
    elseif k == 500
        plot(x, m,'m','linewidth',1.5)
        hold on
    elseif k == 700
        plot(x, m,'g','linewidth',1.5)
        hold on
    elseif k == 1400
        plot(x, m,'c','linewidth',1.5)
        grid on
        xlabel('Non-dimensional Distance through Nozzle','fontsize',12);
        ylabel('Non-dimensional mass flow rate','fontsize',12);
        title('Non-dimensional mass flow rate Vs Distance through nozzle','color','r','fontsize',14);
        legend({'Time step 50','Time step 100','Time step 150','Time step 200','Time step 500','Time step 700','Time step 1400'});
        
    end
end
end

Code Explanation:

1) The function defines the input for the code structure and provided expected outputs. The all underscore t value are at the throat position.

2) The inputs are defined as per the given in the book and time step is defined with the loop from which lowest absolute value is taken as time step for the whole simulation.

3) In the time loop entire program is written with predictor step and corrector step.

4) After finding out the value for next time step the orginal value are calculated and compared with the numerical results.

5) The plot is created for the variation of mass flow rate for various time step numbers.

The main function code for the Conservative form is given as:

function [v,rho,t,m,mach,pr,v_t,rho_t,t_t,m_t,mach_t,pr_t] = conservative_equation(n,L,c,nt,throat)

%Inputs

%inputs
x = linspace(0,L,n);
dx = x(2)-x(1);
a = 1+2.2*(x-1.5).^2;


% Initial conditions Ref J. Anderson Book
for i = 1:n
    if (x(i)>= 0 && x(i)<=0.5)
        rho(i) = 1;
        t(i) = 1;
    elseif (x(i)>= 0.5 && x(i)<=1.5)
        rho(i) = 1-0.366*(x(i)-0.5);
        t(i) = 1-0.167*(x(i)-0.5);
    elseif (x(i)>= 1.5 && x(i)<=3.5)
        rho(i) = 0.634-0.3879*(x(i)-1.5);
        t(i) = 0.833-0.3507*(x(i)-1.5);
    end
end

v = 0.59./(rho.*a);
gamma = 1.4;

% Solution Vectors

pr = rho.*t;
u1 = rho.*a;
u2 = rho.*a.*v;
u3 = u1.*((t./(gamma-1))+((gamma/2).*(v.^2)));

% Time step value
dt = min(c.*(dx./((t.^0.5)+v)));

% Time loop

for k = 1:nt
    
    u1_old = u1;
    u2_old = u2;
    u3_old = u3;
    
    %Flux vectors F
    f1 = u2;
    f = (((gamma-1)/gamma)*(u3-((gamma/2)*((u2.^2)./u1))));
    f2 = ((u2.^2)./u1)+f;
    f3 = ((gamma*u2.*u3)./u1)-(gamma*(gamma-1)*(u2.^3))./(2*u1.^2);
    
    % Predictor step
    
    for i = 2:n-1
        
        df1dx(i) = (f1(i+1)-f1(i))/dx;
        df2dx(i) = (f2(i+1)-f2(i))/dx;
        df3dx(i) = (f3(i+1)-f3(i))/dx;
        dlogadx(i) = (log(a(i+1))-log(a(i)))/dx;
        
        %j_term(i) = f(i)*dlogadx(i);
        dadx(i) = (a(i+1)-a(i))/dx;
        j2(i) = (1/gamma)*(rho(i)*t(i))*dadx(i);
        
        % Continuity equation
        du1dt_p(i) = -df1dx(i);
        % Momentum equation
        du2dt_p(i) = -df2dx(i)+j2(i);
        %Energy equation
        du3dt_p(i) = -df3dx(i);
        
        % Solution updates
        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
    
    % Primitive variable at each time steps
    rho = u1./a;
    v = u2./u1;
    t = (gamma-1)*((u3./u1)-((gamma*v.^2/2)));
    
    %Pure flux terms for corrector step
    
    f1 = u2;
    T1 = (((gamma-1)/gamma)*(u3-(gamma*u2.^2)./(2*u1)));
    f2 = ((u2.^2)./u1)+T1;
    f3 = ((gamma*u2.*u3)./u1)-(gamma*(gamma-1)*(u2.^3))./(2*u1.^2);
    
    % Corrector step
    for i = 2:n-1
        df1dx(i) = (f1(i)-f1(i-1))/dx;
        df2dx(i) = (f2(i)-f2(i-1))/dx;
        df3dx(i) = (f3(i)-f3(i-1))/dx;
        dlogadx(i) = (log(a(i))-log(a(i-1)))/dx;
        
        dadx1(i) = (a(i)-a(i-1))/dx;
        j3(i) = (1/gamma)*(rho(i)*t(i))*dadx1(i);
        
        % Continuity equation
        du1dt_c(i) = -df1dx(i);
        % Momentum equation
        du2dt_c(i) = -df2dx(i)+j3(i);
        %Energy equation
        du3dt_c(i) = -df3dx(i);
        
    end
    
    %Average term derivative
    
    du1dt = 0.5*(du1dt_p + du1dt_c);
    du2dt = 0.5*(du2dt_p + du2dt_c);
    du3dt = 0.5*(du3dt_p + du3dt_c);
    
    % Final time step value
    
    for s = 2:n-1
        u1(s) = u1_old(s) + du1dt(s)*dt;
        u2(s) = u2_old(s) + du2dt(s)*dt;
        u3(s) = u3_old(s) + du3dt(s)*dt;
    end
    
    % Apply boundary conditions
    
    % Inlet BC
    
    u1(1) = rho(1)*a(1);
    u2(1) = 2*u2(2) - u2(3);
    v(1) = u2(1)/u1(1);
    u3(1) = u1(1)*((t(1)/(gamma-1))+((gamma/2)*v(1)^2));
    
    % Outlet BC
    
    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);
    
    % Calculating original terms
    rho = u1./a;
    v = u2./u1;
    t = (gamma-1)*((u3./u1)-((gamma*v.^2)/2));
    
    % Calculation of solution terms
    
    m = rho.*v.*a;
    mach = v./((t).^0.5);
    pr = rho.*t;
    
    %throat value
    
    rho_t(k) = rho(throat);
    t_t(k) = t(throat);
    pr_t(k) = pr(throat);
    mach_t(k) = mach(throat);
    v_t(k) = v(throat);
    m_t(k) = m(throat);
    
    %Plotting mass flow rate
    
    figure(2);
    if k == 50
        plot(x, m,'b','linewidth',1.5)
        hold on
    elseif k==100
        plot(x, m,'r','linewidth',1.5)
        hold on
    elseif k==150
        plot(x, m,'k','linewidth',1.5)
        hold on
    elseif k==200
        plot(x, m,'y','linewidth',1.5)
        hold on
    elseif k == 500
        plot(x, m,'m','linewidth',1.5)
        hold on
    elseif k == 700
        plot(x, m,'g','linewidth',1.5)
        hold on
    elseif k == 1400
        plot(x, m,'c','linewidth',1.5)
        grid on
        xlabel('Non-dimensional Distance through Nozzle','fontsize',12);
        ylabel('Non-dimensional mass flow rate','fontsize',12);
        title('Non-dimensional mass flow rate Vs Distance through nozzle','color','r','fontsize',14);
        legend({'Time step 50','Time step 100','Time step 150','Time step 200','Time step 500','Time step 700','Time step 1400'});
    end
end
end

Code Explanation:

1) For the conservative form some different approach is used to write the differential equation to find out the flow variables.

2) To ease the code and lower the amount of debugging in the code we are defining some terms which will be assosiated with the general equation. Defined in the tems of solution vectors as 'U' and flux vectors 'f' with source term 'j'.

3) Then this distributed form are used to write the complete code developement for each equations. Then similar to non-conservative form the predictor step and corrector steps are being calculated to get the next time step value.

4) After that the throat value are calculated and plot is made to determine the mass flow rate at different time steps value.

Results:

1) Mass flow rate through the Nozzle distance for Non-conservative form:

A) As we can see from the graph for time step number after 500 the solution seems to be more stable and overlaping each other.

B) Hence we can see that throughout the nozzle distance the solution get steady and dont have any fluctuations. The mass flow rate at which we are getting the the stable solution is 0.6.

2) Variation of flow parameters along x axis

As we have stated earlier we are simulating the 1D supersonic nozzle flow, hence it is important to capture the parameters which changes throughout that dimension.

As we can see in the above graph the value for each primitive variable are approximately matching with exact solution for both conservative and non-conservative forms of equation.

3) Time variation of flow parameters

A) The above plot shows us the timewise variation of primitive parameters through the solution time. If we compare the acutal solution for the both conservative and non-conservative form the it is matches exactly.

B) From the graph it is concluded that the non-conservative solution provide the approximation solution near to the exact solution.

C) Also we can see that the curve gains the steady state position after 500 time steps for both conservative and non-coservative form of equations.

4) Mass flow rate through the Nozzle distance for Conservative form:

A) The fluctuation in above plot are less as compare to the plot of curve for Non-conservative part.

B) The fluctuation are not that much important as we have selected different initial condition for the conservative form equation with constant mass flow rate.

C) The same behaviour is seen for the time step number 500,700 and 1400 as they are overlapping the same region of space compared to the non-conservative part of graph.

D) Compared to the non-conservative form the conservative form of equation gives the better results when compare with analytical solution which can be seen from the below graph.

E) From the above graph we can see clearly that conservative form steadier and more accurate as compare to the non-conservative form.

F) The only reason is that, in the conservative form mass flow rate is one of the dependent variable comes in the equation while in non-conservative form the mass flow rate is calculated from the dependent variables or primitive variables.

5) Error calculation for Density, Temperature, Velocity and Mass Flow rate

For the calculation of error in the non-conservative form and conservative form we will introduce the new code but just with some modification. As it was not included in the original main code because it was creating limited number of time step number because of that for each value of tolerance it need to be change accordingly which creating more debugging of code structure.

The main Code is given as:

clear all
close all
clc

n = 31;                     % No. of gris points
L = 3;                      % Length of Domain
x = linspace(0,L,n);        % Discretization
c = 0.5;                    % Value of courant number
throat = ((n-1)/2)+1;       % Throat dimension
nt = 1400;                  % Number of time steps
tol = 1e-6;


% Non-conservative form of Equation
tic
[v,rho,t,m1,mach,pr,v_t,rho_t,t_t,m_t,mach_t,pr_t,er_v1,er_rho1,er_t1,er_mass1] = non_conservative_equation1(n,c,L,nt,throat,tol);
Non_Conservative_time = toc;

tic
[v2,rho2,t2,m2,mach2,pr2,v_t2,rho_t2,t_t2,m_t2,mach_t2,pr_t2,er_v2,er_rho2,er_t2,er_mass2] = conservative_equation2(n,L,c,nt,throat,tol);
Conservative_time = toc;

%Error plots
figure(6)
subplot(2,1,1)
plot(er_v1);
axis([0 inf -0.01 0.04]);
xlabel('time step')
ylabel('Error in velocity')
title('Error in velocity for Non-coservative form')
grid on
subplot(2,1,2)
plot(er_v2)
axis([0 inf -0.01 0.04]);
xlabel('Time step')
ylabel('Error in velocity')
title('Error in velocity for Conservative form')
grid on

figure(7)
subplot(2,1,1)
plot(er_t1);
axis([0 inf -0.005 0.01]);
xlabel('time step')
ylabel('Error in Temperature')
title('Error in Temperature for Non-coservative form')
grid on
subplot(2,1,2)
plot(er_t2)
axis([0 inf -0.005 0.01]);
xlabel('Time step')
ylabel('Error in Temperature')
title('Error in Temperature for Conservative form')
grid on

figure(8)
subplot(2,1,1)
plot(er_rho1);
axis([0 inf -0.01 0.02]);
xlabel('time step')
ylabel('Error in Density')
title('Error in Density for Non-coservative form')
grid on
subplot(2,1,2)
plot(er_rho2)
axis([0 inf -0.01 0.02]);
xlabel('Time step')
ylabel('Error in Density')
title('Error in Density for Conservative form')
grid on

figure(9)
subplot(2,1,1)
plot(er_mass1);
xlabel('time step')
ylabel('Error in Mass flow rate')
title('Error in Mass flow rate for Non-coservative form')
grid on
subplot(2,1,2)
plot(er_mass2)
xlabel('Time step')
ylabel('Error in Mass flow rate')
title('Error in Mass flow rate for Coservative form')
grid on

The function for Non-conservative form is given as:

function [v,rho,t,m1,mach,pr,v_t,rho_t,t_t,m1_t,mach_t,pr_t,er_rho,er_v,er_t,er_mass] = non_conservative_equation1(n,c,L,nt,throat,tol)

%inputs
x = linspace(0,L,n);
dx = x(2)-x(1);

%Initial conditions
rho = 1-0.3146*x;
t = 1-0.2314*x;
v = (0.1+1.09*x).*t.^0.5;
a = 1+2.2*(x-1.5).^2;
gamma = 1.4;

%Time step value
for i = 1:n
    b(i) = c*(dx/(t(i)^0.5+v(i)));
end
dt = min(abs(b));

% Time Loop
for k = 1:nt
    
    rho_old = rho;
    v_old = v;
    t_old = t;
    
    %Predictor methods
    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; %keep track of this bracket
        
        %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 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; %keep track of this bracket
        
        %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
    %Average derivative
    drhodt = 0.5*(drhodt_p + drhodt_c);
    dvdt = 0.5*(dvdt_p + dvdt_c);
    dtdt = 0.5*(dtdt_p + dtdt_c);
    
    % Final time step
    for p = 2:n-1
        rho(p)= rho_old(p) + drhodt(p)*dt;
        v(p) = v_old(p) + dvdt(p)*dt;
        t(p) = t_old(p) + dtdt(p)*dt;
    end
    
    % Apply boundary condtions
    % Inlet BC
    v(1) = 2*v(2)-v(3);
    rho(1) = 1; % Check this
    t(1) = 1; % Check this also
    
    %Outlet BC
    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);
    
    % Calculation of solutions
    
    m1 = rho.*v.*a;
    mach = v./((t).^0.5);
    pr = rho.*t;
    
    %throat value
    
    rho_t(k) = rho(throat);
    t_t(k) = t(throat);
    pr_t(k) = pr(throat);
    mach_t(k) = mach(throat);
    v_t(k) = v(throat);
    m1_t(k) = m1(throat);
    
    % Error in mass flow rate
    inflow = rho(1)*v(1)*a(1);
    outflow = rho(n)*v(n)*a(n);
    mass_netflow = inflow - outflow;
    
    % Error in the primitive value 
    er_v(k) = max(v-v_old);
    er_rho(k) = max(rho-rho_old);
    er_t(k) = max(t - t_old);
    er_mass(k) = max(m1 - mass_netflow);
    if (er_v(k)

The function for Conservative form is given as:

function [v,rho,t,m,mach,pr,v_t,rho_t,t_t,m_t,mach_t,pr_t,er_v,er_rho,er_t,er_mass] = conservative_equation2(n,L,c,nt,throat,tol)

%Inputs

%inputs
x = linspace(0,L,n);
dx = x(2)-x(1);
a = 1+2.2*(x-1.5).^2;


% Initial conditions Ref J. Anderson Book
for i = 1:n
    if (x(i)>= 0 && x(i)<=0.5)
        rho(i) = 1;
        t(i) = 1;
    elseif (x(i)>= 0.5 && x(i)<=1.5)
        rho(i) = 1-0.366*(x(i)-0.5);
        t(i) = 1-0.167*(x(i)-0.5);
    elseif (x(i)>= 1.5 && x(i)<=3.5)
        rho(i) = 0.634-0.3879*(x(i)-1.5);
        t(i) = 0.833-0.3507*(x(i)-1.5);
    end
end

v = 0.59./(rho.*a);
gamma = 1.4;

% Solution Vectors

pr = rho.*t;
u1 = rho.*a;
u2 = rho.*a.*v;
u3 = u1.*((t./(gamma-1))+((gamma/2).*(v.^2)));

% Time step value
dt = min(c.*(dx./((t.^0.5)+v)));

% Time loop

for k = 1:nt
    
    u1_old = u1;
    u2_old = u2;
    u3_old = u3;
    
    %Flux vectors F
    f1 = u2;
    f = (((gamma-1)/gamma)*(u3-((gamma/2)*((u2.^2)./u1))));
    f2 = ((u2.^2)./u1)+f;
    f3 = ((gamma*u2.*u3)./u1)-(gamma*(gamma-1)*(u2.^3))./(2*u1.^2);
    
    % Predictor step
    
    for i = 2:n-1
        
        df1dx(i) = (f1(i+1)-f1(i))/dx;
        df2dx(i) = (f2(i+1)-f2(i))/dx;
        df3dx(i) = (f3(i+1)-f3(i))/dx;
        dlogadx(i) = (log(a(i+1))-log(a(i)))/dx;
        
        %j_term(i) = f(i)*dlogadx(i);
        dadx(i) = (a(i+1)-a(i))/dx;
        j2(i) = (1/gamma)*(rho(i)*t(i))*dadx(i);
        
        % Continuity equation
        du1dt_p(i) = -df1dx(i);
        % Momentum equation
        du2dt_p(i) = -df2dx(i)+j2(i);
        %Energy equation
        du3dt_p(i) = -df3dx(i);
        
        % Solution updates
        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
    
    % Primitive variable at each time steps
    rho = u1./a;
    v = u2./u1;
    t = (gamma-1)*((u3./u1)-((gamma*v.^2/2)));
    
    %Pure flux terms for corrector step
    
    f1 = u2;
    T1 = (((gamma-1)/gamma)*(u3-(gamma*u2.^2)./(2*u1)));
    f2 = ((u2.^2)./u1)+T1;
    f3 = ((gamma*u2.*u3)./u1)-(gamma*(gamma-1)*(u2.^3))./(2*u1.^2);
    
    % Corrector step
    for i = 2:n-1
        df1dx(i) = (f1(i)-f1(i-1))/dx;
        df2dx(i) = (f2(i)-f2(i-1))/dx;
        df3dx(i) = (f3(i)-f3(i-1))/dx;
        dlogadx(i) = (log(a(i))-log(a(i-1)))/dx;
        
        dadx1(i) = (a(i)-a(i-1))/dx;
        j3(i) = (1/gamma)*(rho(i)*t(i))*dadx1(i);
        
        % Continuity equation
        du1dt_c(i) = -df1dx(i);
        % Momentum equation
        du2dt_c(i) = -df2dx(i)+j3(i);
        %Energy equation
        du3dt_c(i) = -df3dx(i);
        
    end
    
    %Average term derivative
    
    du1dt = 0.5*(du1dt_p + du1dt_c);
    du2dt = 0.5*(du2dt_p + du2dt_c);
    du3dt = 0.5*(du3dt_p + du3dt_c);
    
    % Final time step value
    
    for s = 2:n-1
        u1(s) = u1_old(s) + du1dt(s)*dt;
        u2(s) = u2_old(s) + du2dt(s)*dt;
        u3(s) = u3_old(s) + du3dt(s)*dt;
    end
    
    % Apply boundary conditions
    
    % Inlet BC
    
    u1(1) = rho(1)*a(1);
    u2(1) = 2*u2(2) - u2(3);
    v(1) = u2(1)/u1(1);
    u3(1) = u1(1)*((t(1)/(gamma-1))+((gamma/2)*v(1)^2));
    
    % Outlet BC
    
    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);
    
    % Calculating original terms
    rho = u1./a;
    v = u2./u1;
    t = (gamma-1)*((u3./u1)-((gamma*v.^2)/2));
    
    % Calculation of solution terms
    
    m = rho.*v.*a;
    mach = v./((t).^0.5);
    pr = rho.*t;
    
    %throat value
    
    rho_t(k) = rho(throat);
    t_t(k) = t(throat);
    pr_t(k) = pr(throat);
    mach_t(k) = mach(throat);
    v_t(k) = v(throat);
    m_t(k) = m(throat);
    
    % Error in mass flow rate
    inflow = rho(1)*v(1)*a(1);
    outflow = rho(n)*v(n)*a(n);
    mass_netflow = inflow - outflow;
    
    %Plotting mass flow rate
    
    % Error in the primitive value 
    er_v(k) = max(u2-u2_old);
    er_rho(k) = max(u1-u1_old);
    er_t(k) = max(u3 - u3_old);
    er_mass(k) = max(m - mass_netflow);
    if (er_v(k)

In both of function the changes is calculation of error of primitive values which is shown accordingly in the end of both functions.

The results are plotted were given as:

A) We can see from the above four graphs for non-conservative form the error gets stabilize at earlier time steps rather than the conservative form.

B) If we lower the tolerance value again the stabilization for the non-conservative form will be higher as compare to the conservative form.

C) In case of mass flow rate for conservative form it shows the lower error from start of the time step whereas in the non-conservative form it is not stabilize and for specific time step it has more error percentage.

6) Check for grid dependency:

The grid independency is one of the most important aspect in the all kind of CAE simulation. It mostly keen toward to the accuracy of the solution and it stability by changing the grid size or the time step size. As we are dealing with 1D simulation the grid dependency check is not require that much if you have made the correct assumption and choose correct boundary conditions, but still it is a great practice by an engineer to make it all side clear. The grid dependency is necessary with 3D elements as the normal of surface changes abruptly depending upon the body shape and size.

The grid dependency can be done by changing the value of grid size and checking how much the primitive variable or dependent variable are changing.

In our case we will do the changes for supersonic 1D nozzle simulation by changing the grid size to 61 which is double the original one for which whole code was written.

We have to also take care for the time step size as it is also depend on the CFL number we have take for each grid points and then we have choose the min value of that as the 'dt'.

The time step size for n=31 is 0.0201 for non-conservative form and for conservative form it is 0.02036. So if we calculate the time elasped for the non-conservative form it is 0.0201*1400 = 28.14 and similarly for conservative form it is 28.504. Hence if we take the n=61 the dt becomes for conservative form as 0.01018 and the time step size will be equal to 2800.

We will compare the value at the throat for n=31 and n=61. Also we will check the respective graph for those conditions.

Parameters	Non-Conservative form		Conservative form		Exact Result
	N = 31	N = 61	N = 31	N = 61
rho_t	0.6387	0.6374	0.6524	0.6391	0.634
pr_t	0.5342	0.5324	0.5487	0.5340	0.528
t_t	0.8365	0.8353	0.8410	0.8356	0.833
mach_t	0.9994	0.9998	0.9780	0.9978	1.00

As we can see the value for the n=61 are nearly matching with the n=31 data set hence we can say that n=31 is grid independent.

We will also check the graph which will generate for the n=61.

A) Timewise variation of parameter for n= 61

B) Variation of Flow parameter

Hence, from above two graph it is clearly seen that there is no significant change in the timewise variation of flow parameter at throat and also variation of flow parameter along x axis.

C) Mass flow rate at different time step value for Non-conservative form at n=61.

As we can see from the above figure and compared with n=31 non-conservative mass flow rate graph increasing in number of grid point creates the better results at earlier time step. In above graph we are getting the steady state solution at time step 1000 when compared with the n=31 non-conservative form graph it was at 500. Hence we can say that n=31 grid size is independent.

D) Mass flow rate at different time step for Conservative form at n=61

E) Comparison of Mass flow rate for Conservative and Non-conservative form at n=61.

Compared with mass flow rate of n=31 we can see that for non-conservative form the error is decrease is marginally and become steady state for throat.

7) Time taken for the simulation

The time required to calculate the Non-conservative form with n=31 is 2.2898 seconds.

The time required to calculate the conservative form with n=31 is 3.0267 seconds.

From the time required to complete the simulation we can say that Non-conservative form is faster than conservative form.

The reason behind the conservative form is slow because of flux variable are used to calculate the primitive variable where as it is not required in case of non-conservative form.