Non-cons vs Cons

A. Main script

% Main program 
clear all
close all
clc

% inputs variables
n = 61;
L = 3;
x = linspace(0,L,n);

% Courant number
c = 0.5;

% Throat Node
throat = ((n-1)/2)+1;

% Time steps
nt = 1400;


% Non conservative form
[m,computational_time_nc,ind,grid_ind_param] = non_conservative_eq(n,L,c,nt,throat);

% Conservative form
[m_c,computational_time_c,ind_c,grid_ind_param_c] = conservative_eq(n,L,c,nt,throat);

% Mass flow rate comparison
figure(7)
plot(x,m,\'r\')
hold on
plot(x,m_c,\'b\')
hold on
plot(x,linspace(0.579,0.579,n),\'--k\')
ylabel(\'Non-dimensional mass flow rate\');
xlabel(\'Non-dimensional distance through nozzle\');
legend({\'Non-Conservative Form\', \'Conservative Form\',\'Exact analytical solution\'});
title({\'Non-dimensional mass flow rate\',\'Non-Conservative Form vs Conservative Form\'}); 
axis([0 3 0.575 0.600]); grid on

B. Non_conservative_eq

function [m,computational_time,ind,grid_ind_param] = non_conservative_eq(n,L,c,nt,throat)
% [m,computational_time] = non_conservative_flow(n,L,c,nt,throat)
tic

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

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

% Time step based on CFL condition
dt = min(c*(dx./(sqrt(t)+v)));

% Outer time loop
for time = 1:nt
    % Saving old values
    rho_old = rho;
    v_old = v;
    t_old = t;

    % 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;
    
        drhodt_p(j) = -rho(j)*dvdx-rho(j)*v(j)*dlogadx-v(j)*drhodx;    % continuity equation
        dvdt_p(j) = -v(j)*dvdx-(1/gamma)*(dtdx +(t(j)/rho(j))*drhodx); % momentum equation
        dtdt_p(j) = -v(j)*dtdx-(gamma-1)*t(j)*(dvdx+v(j)*dlogadx);     % energy equation

        % 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;

        drhodt_c(j) = -rho(j)*dvdx -rho(j)*v(j)*dlogadx -v(j)*drhodx;   % continuity equation
        dvdt_c(j) = -v(j)*dvdx -(1/gamma)*(dtdx +(t(j)/rho(j))*drhodx); % momentum equation
        dtdt_c(j) = -v(j)*dtdx - (gamma-1)*t(j)*(dvdx + v(j)*dlogadx);  % energy equation
    end 
    
    % Compute 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 j = 2: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);
    % rho(1) and t(1) are fixed and independent of time
  
    % 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);
  
    % Calculating other variables
    m = rho.*v.*a;        % Mass flow rate
    mach = v./((t).^0.5); % Mach number
    pr = rho.*t;          % Pressure ratio
  
    % Time series at throat at time = 1400dt
    rho_t(time) = rho(throat);
    t_t(time) = t(throat);
    pr_t(time) = pr(throat);
    mach_t(time) = mach(throat);
    v_t(time) = v(throat);
    m_t(time) = m(throat);
    
    % Plot of mass flow rate at different time steps
    figure(1);
    if time  == 1
        plot(x,m,\'y\')
        hold on
    elseif time == 50 
        plot(x,m,\'m\')
        hold on
    elseif time == 100
        plot(x,m,\'c\')
        hold on
    elseif time == 150
        plot(x,m,\'r\')
        hold on
    elseif time == 200
        plot(x,m,\'g\')
        hold on
    elseif time == 550
        plot(x,m,\'b\') 
        hold on 
    elseif time == 700
        plot(x,m,\'w\')
        hold on
    elseif time == 1400
        plot(x,m,\'k\')
        grid on
        ylabel(\'Non-dimensional mass flow rate\');
        xlabel(\'Non-dimensional distance through nozzle\');
        legend({\'1dt\',\'50dt\',\'100dt\',\'150dt\',\'200dt\',\'550dt\',\'700dt\',\'1400dt\'});
    end  
end

computational_time = toc;

% Grid dependance
grid_ind_param = [rho_t(time),t_t(time),pr_t(time),mach_t(time),v_t(time),m_t(time)];

% Min time step to converge
for ii = 1:nt-1
    error(ii) = abs(rho_t(ii+1)-rho_t(ii))/rho_t(ii+1)*100;
    if error(ii) < 0.001
        ind = ii;
        break
    end
end

% Plots
figure(2)
plot1 = subplot(4,1,1);
plot(linspace(1,nt,nt),rho_t); grid on
title({\'Timewise variation of non-dimensional Density\',\'Temperature , Pressure and Mach number at throat\'});
plot2 = subplot(4,1,2);
plot(linspace(1,nt,nt),t_t); grid on
plot3 = subplot(4,1,3);
plot(linspace(1,nt,nt),pr_t); grid on
plot4 = subplot(4,1,4);
plot(linspace(1,nt,nt),mach_t); grid on

figure(3)
plot1 = subplot(4,1,1);
plot(linspace(1,n,n),rho);grid on; 
title({\'Steady-state distributions of non-dimensional Density\',\'Temperature , Pressure and Mach number\',\'as function of non-dimensional distance\'});
plot2 = subplot(4,1,2);
plot(linspace(1,n,n),t); grid on
plot3 = subplot(4,1,3);
plot(linspace(1,n,n),pr); grid on
plot4 = subplot(4,1,4);
plot(linspace(1,n,n),mach); grid on

end

C. Conservative_eq

function [m,computational_time,ind,grid_ind_param] = conservative_eq(n,L,c,nt,throat)
% [m,computational_time] = conservative_flow(n,L,c,nt,throat)
tic

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

% Calculate improved initial profiles
for i = 1:n
    if (x(i)>=0 && x(i)<=0.5)
        rho(i) = 1; % Density
        t(i) = 1;   % Temperature
    elseif (x(i)>=0.5 && x(i)<=1.5)
        rho(i) = 1-0.366*(x(i)-0.5); % Density
        t(i) = 1-0.167*(x(i)-0.5);   % Temperature
    elseif (x(i)>=1.5 && x(i)<=3.5)
        rho(i) = 0.634-0.3879*(x(i)-1.5); % Density
        t(i) = 0.833-0.3507*(x(i)-1.5);   % Temperature
    end
end
a = 1+2.2*(x - 1.5).^2; % Area 
v = 0.59./(rho.*a);     % Velocity (Assuming steady state flow)
  
% Solution Vector U
U1 = rho.*a;
U2 = rho.*a.*v;
U3 = U1.*((t./(gamma-1))+((gamma/2).*(v.^2)));

% Time step based on CFL condition
dt = min(c.*(dx./((t.^0.5)+v))); 
 
% Outer time loop
for time = 1:nt
    % Saving old values  
    U1_old = U1;
    U2_old = U2;
    U3_old = U3;
    
    % PURE FORM OF THE FLUX TERMS
    F1 = U2;
    term = (((gamma-1)/gamma)*(U3-((gamma/2)*((U2.^2)./U1))));
    F2 = ((U2.^2)./U1)+ term;
    F3 = ((gamma*U2.*U3)./U1)-(gamma*(gamma-1)*(U2.^3))./(2*U1.^2);

    % Predictor Method
    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;
        dadx(i) = (a(i+1)-a(i))/dx;
        J2(i) = (1/gamma)*(rho(i)*t(i))*dadx(i); 
        
        dU1dt_p(i) = -dF1dx(i);         % Continuity equation
        dU2dt_p(i) = -dF2dx(i) + J2(i); % Momentum equation
        dU3dt_p(i) = -dF3dx(i);         % Energy equation
        
        % 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
    
    % Initialize the values for the corrector step
    rho = U1./a;
    v = U2./U1;
    t = (gamma-1)*((U3./U1)-((gamma*v.^2)/2));
    
    % PURE FORM OF THE FLUX TERMS
    F1 = U2;
    term1 = (((gamma-1)/gamma)*(U3-(gamma*U2.^2)./(2*U1)));
    F2 = ((U2.^2)./U1)+ term1;
    F3 = ((gamma*U2.*U3)./U1)-(gamma*(gamma-1)*(U2.^3))./(2*U1.^2);
    
    % Corrector Method
    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); 
 
        
        dU1dt_c(i) = -dF1dx(i);         % Continuity equation
        dU2dt_c(i) = -dF2dx(i) + J3(i); % Momentum equation
        dU3dt_c(i) = -dF3dx(i);         % Energy equation  
    end
    
    % Average derivative
    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: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
    
    % Boundary conditions
    % Inlet
    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 
    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);
    
    % Final solution update
    rho = U1./a; 
    v = U2./U1;
    t = (gamma-1)*((U3./U1)-((gamma*v.^2)/2));
    
   
    % Calculating other variables 
    m = rho.*v.*a;
    mach = v./((t).^0.5);
    pr = rho.*t;
    
    % Throat values 
    rho_t(time) = rho(throat);
    t_t(time) = t(throat); 
    pr_t(time) = pr(throat);
    mach_t(time) = mach(throat);
    v_t(time) = v(throat);
    m_t(time) = m(throat);
    
    % Plot of mass flow rate at different time steps
    figure(4);
    if time == 1
        plot(x,m,\'y\')
        hold on
    elseif time == 100
        plot(x,m,\'m\')
        hold on
    elseif time == 200
        plot(x,m,\'c\')
        hold on
    elseif time == 700
        plot(x,m,\'r\')
        hold on
    elseif time  == 1400
        plot(x,m,\'g\')
        grid on
        ylabel(\'Non-dimensional mass flow rate\');
        xlabel(\'Non-dimensional distance through nozzle\');
        legend({\'1dt\',\'100dt\',\'200dt\',\'700dt\',\'1400dt\'});
    end
end

computational_time = toc;

% Grid dependance
grid_ind_param = [rho_t(time),t_t(time),pr_t(time),mach_t(time),v_t(time),m_t(time)];

% Min time step to converge
for ii = 1:nt-1
    error(ii) = abs(rho_t(ii+1)-rho_t(ii))/rho_t(ii+1)*100;
    if error(ii) < 0.001
        ind = ii;
        break
    end
end

% Plots
figure(5)
plot1 = subplot(4,1,1);
plot(linspace(1,nt,nt),rho_t); grid on
title({\'Timewise variation of non-dimensional Density\',\'Temperature , Pressure and Mach number at throat\'});
plot2 = subplot(4,1,2);
plot(linspace(1,nt,nt),t_t); grid on
plot3 = subplot(4,1,3);
plot(linspace(1,nt,nt),pr_t); grid on
plot4 = subplot(4,1,4);
plot(linspace(1,nt,nt),mach_t); grid on

figure(6)
plot1 = subplot(4,1,1);
plot(linspace(1,n,n),rho);grid on; 
title({\'Steady-state distributions of non-dimensional Density\',\'Temperature , Pressure and Mach number\',\'as function of non-dimensional distance\'});
plot2 = subplot(4,1,2);
plot(linspace(1,n,n),t); grid on
plot3 = subplot(4,1,3);
plot(linspace(1,n,n),pr); grid on
plot4 = subplot(4,1,4);
plot(linspace(1,n,n),mach); grid on
end

With regards to this algorithm, there are three different parts which satisfy the requirements of the challenge. The first part of the script is the main programme in which necessary parameters are defined, two different functions are called, and a figure of non-dimensional mass flow rate comparison is plotted. Then, the non-conservative equations are solved in a single function. Its outputs are the mass flow rate, the computational time, the min number of iterations for the error to be less than 0.001%, a matrix with some flow parameters (for the grid dependence test) and three plots: 1) The Non-dimensional mass flow rate for different time steps, 2) The timewise variation of non-dimensional flow parameters at throat and 3) The steady-state distributions of non-dimensional flow parameter as function of non-dimensional distance. Same outputs for the function which solves the conservative equations.

With regards to the min iterations for convergence, Non-conservative equations need 376 iterations while conservative 247 iterations. The total computational time for non-conservation equations to be solved is 4.36sec while for the conservative eq is 4.1 sec. The following table shows the values of some flow parameters for 31 points and 61 points of the computational grid.

	Density	Temperature	Pressure	Mach Number	Velocity	Mass flow
31p Non-cons	0,638664	0,83645072	0,534211	0,999387926	0,914017	0,58375
61p Non-cons	0,637659	0,835439136	0,532726	0,999775569	0,913818	0,5827049
31p Cons	0,652367	0,841029479	0,54866	0,978016215	0,896916	0,5851186
61p Cons	0,639077	0,835549624	0,533981	0,997851322	0,91212	0,582915

It should be mentioned that all the algorithm is based on the 7th chapter of the “CFD, the basics with applications” Jonh D Anderson