Challenge

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

Aim:

  • To simulate the isentropic flow through a quasi 1D subsonic-supersonic nozzle.
  • To derive both the conservation and non-conservation forms of the governing equations and sovle them using the MacCormack's technique.

Problem statement :

  • To determine the steady-state temperature distribution for the flow-field variables and to investigate the difference between the two forms of governing equations by comparing their solutions.
  • To plot steady-state distribution of primitive variables inside the nozzle.
  • To plot time-wise variation of the primitive variables.
  • To plot variation of Mass flow rate distribution inside the nozzle at different time steps during the time-marching process.
  • To plot comparison of Normalized mass flow rate distributions of both forms.

Nozzle diagram :

  • At convergent section the flow velocity is less than the sound velocity therefore there mach number is less than 1.
  • At throat section the flow velocity is equal to the sound velocity therefore there mach number is equal to 1.
  • At divergent section the flow velocity is greater than the sound velocity therefore there mach number is greater than 1.

Non Dimensional Quantities :

Also as part of the problem setup we represent quantities in terms of non-dimensional quantities

Non Conservative form :

The non-conservative form has its differentiated variables as coefficients of some derivatives in the equation.

Boundary conditions :

In our case, we are having a subsonic flow in the convergent section (M<1), choked flow at the throat section (M=1) and supersonic outflow in the divergent section (M>1). So by using method of characteristics for an unsteady, one-dimensional flow we can precisely allot inflow and outflow boundary conditions and can identify which quantities to state at the boundaries.

Conservative form :

In conservative form, we will initialize the nozzle profile rationally. As in conservative form, seen above from the governing equations, we have solution vectors, flux terms and a source term. So it is the prior matter of concern to state initial and boundary conditions in those terms rather than that of primitive variables.

Governing Equations
Generic Form of Governing Equations
Initial Conditions
Boundary Conditions:
CODE FOR NON-CONSERVATIVE FORM
clc 
clear all
close all

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

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

% Time Step
nt = 1400;
CFL = 0.5;

th = round(n/2);

% Time Loop
for i = 1:nt
     
    dt(i) = min((CFL*dx)./(v+sqrt(T)));
    
    rho_old = rho;
    T_old = T;
    v_old = v;
       
    % Predictive Method
    for j = 2:n-1
        
        dvdx = (v(j+1)-v(j))/dx;
        dlogA_dx = (log(A(j+1))-log(A(j)))/dx;
        drhodx  = (rho(j+1)-rho(j))/dx;
        dTdx = (T(j+1)-T(j))/dx;
        
        %Prdictive step
        
        % Continuity Equation
        drhodt_p(j) = -rho(j)*dvdx-rho(j)*v(j)*dlogA_dx-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)*dlogA_dx);
     
        %Predictive Solution update
        rho(j) = rho_old(j) + drhodt_p(j)*dt(i);
        v(j) = v_old(j) + dvdt_p(j)*dt(i);
        T(j) = T_old(j) + dTdt_p(j)*dt(i);
    end
    
    % Corrective Method
    
    for j = 2:n-1
        
        dvdx = (v(j)-v(j-1))/dx;
        dlogA_dx = (log(A(j))-log(A(j-1)))/dx;
        drhodx  = (rho(j)-rho(j-1))/dx;
        dTdx = (T(j)-T(j-1))/dx;
        
        %Prdictive step
        
        % Continuity Equation
        drhodt_c(j) = -rho(j)*dvdx-rho(j)*v(j)*dlogA_dx-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)*dlogA_dx);
     
        end
   
    % Average Values
    
    drhodt = 0.5*(drhodt_p+drhodt_c);
    dvdt = 0.5*(dvdt_p+dvdt_c);
    dTdt = 0.5*(dTdt_p+dTdt_c);
    
    % Final Solution
    for k = 2:n-1
        
        rho(k) = rho_old(k) + drhodt(k)*dt(i);
        v(k) = v_old(k) + dvdt(k)*dt(i);
        T(k) = T_old(k) + dTdt(k)*dt(i);    
    end
    
    % BCs
    % Inlet
    
    rho(1) = 1;
    T(1) = 1;
    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);
    
    % Pressure
    p = rho.*T;
    
    % Mach Number 
    M = v./sqrt(T);
    
    % Mass Flow Rate
    m = rho.*A.*v;
    
    % Primitive Variables at throat 
    M_throat(i) = M(th);
    rho_throat(i) = rho(th);
    v_throat(i) = v(th);
    T_throat(i) = T(th);
    p_throat(i) = p(th);
    
    % Mass flow rate at different time steps
    figure(4)
    hold on
    if i ==1
        plot(x,m)
    end
    if i ==50
        plot(x,m)
    end
    if i ==100
        plot(x,m)
    end
    if i ==200
        plot(x,m)
    end
    if i ==400
        plot(x,m)
    end
    if i ==800
        plot(x,m)
    end
             
end

% Plotting 

% Mass flow rate at different time steps
xlabel('Nozzle Length')
ylabel('Mass Flow Rate')
legend('1 dt','50 dt','100 dt','200 dt','400 dt','800 dt')
title('Variation of Mass Flow Rate at different time steps')
% Steady State Primitive Variable Profiles

figure(1)
subplot(5,1,1)
plot(x,M,'b')
xlabel('Nozzle Length')
ylabel('Mach Number')
title('Variation of Steady State Primitive Variable Profiles')
subplot(5,1,2)
plot(x,rho,'k')
xlabel('Nozzle Length')
ylabel('Density')
subplot(5,1,3)
plot(x,v,'m')
xlabel('Nozzle Length')
ylabel('velocity')
subplot(5,1,4)
plot(x,T,'r')
xlabel('Nozzle Length')
ylabel('Temperature')
subplot(5,1,5)
plot(x,p,'y')
xlabel('Nozzle Length')
ylabel('Pressure')

% Transient State Primitive Variable Profiles
figure(2)
subplot(5,1,1)
plot(M_throat)
xlabel('Time Step')
ylabel('Mach Number')
title('Variation of Transient State Primitive Variable Profiles at Nozzle Throat')
subplot(5,1,2)
plot(rho_throat,'k')
xlabel('Time Step')
ylabel('Density')
subplot(5,1,3)
plot(v_throat,'r')
xlabel('Time Step')
ylabel('velocity')
subplot(5,1,4)
plot(T_throat,'g')
xlabel('Time Step')
ylabel('Temperature')
subplot(5,1,5)
plot(p_throat,'y')
xlabel('Time Step')
ylabel('Pressure')

% Mass Flow rate
figure(3)
plot(x,m)
xlabel('Nozzle Length')
ylabel('Mass Flow Rate')
title('Mass Flow Rate along Nozzle length')
grid on

Output of Non Conservative:
MATLAB CODE FOR CONSERVATIVE FORM:
clc 
clear all
close all

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

th = round(n/2);

for i = 1:length(x)
if x(i) <= 0.5 
    rho(i) = 1;
    T(i) = 1;
end
end

for i = 1:length(x)
if 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);
end
end

for i = 1:length(x)
if x(i) >= 1.5
    rho(i) = 0.634-0.3879*(x(i)-1.5);
    T(i) = 0.833-0.3507*(x(i)-1.5);
end
end

% Initial Profiles
A = 1+ 2.2*(x-1.5).^2;
v = 0.59./(rho.*A);

gamma = 1.4;

% Time Step
nt = 1400;
CFL = 0.5;
 
dt = min((CFL*dx)./(v+sqrt(T)));

U1=rho.*A;
U2=rho.*A.*v;
U3=rho.*((T./(gamma-1))+(gamma/2).*v.^2).*A;



% Time Loop
for i = 1:nt
     
       
    U1_old = U1;
    U2_old = U2;
    U3_old = U3;
    
    F1=U2;
    F2=((U2.^2)./U1)+((gamma-1)/gamma)*(U3-(gamma/2)*(((U2).^2)./U1));
    F3=(((gamma)*(U2.*U3))./U1)-(((gamma*(gamma-1))/2)*((U2.^3)./U1.^2));
       
    % Predictive Method
    for j = 2:n-1
        
             
       dU1dt_p(j) = -(F1(j+1)-F1(j))/dx;
        
       J2 = (1/gamma)*rho(j)*T(j)*(A(j+1)-A(j))/dx;
       
       dU2dt_p(j) = -(F2(j+1)-F2(j))/dx + J2;
       
       dU3dt_p(j) = -(F3(j+1)-F3(j))/dx;
           
        %Predictive 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
    
    %Update F & Pritimitive variables vector for J2
    rho = U1./A;
    T = (gamma-1)*(U3./U1-0.5*gamma*(U2./U1).^2);
    
    F1=U2;
    F2=((U2.^2)./U1)+((gamma-1)/gamma)*(U3-(gamma/2)*(((U2).^2)./U1));
    F3=(((gamma)*(U2.*U3))./U1)-(((gamma*(gamma-1))/2)*((U2.^3)./U1.^2));
    
% Corrective Method
    
    
    for j = 2:n-1        
             
             
       dU1dt_c(j) = -(F1(j)-F1(j-1))/dx;
        
       J2 = (1/gamma)*rho(j)*T(j)*(A(j)-A(j-1))/dx;
      
       dU2dt_c(j) = -(F2(j)-F2(j-1))/dx+J2;
       
       dU3dt_c(j) = -(F3(j)-F3(j-1))/dx;
        
            
    end
   
    % Average Values
    
    dU1dt = 0.5*( dU1dt_p+ dU1dt_c);
    dU2dt = 0.5*( dU2dt_p+ dU2dt_c);
    dU3dt = 0.5*( dU3dt_p+ dU3dt_c);
    
    % Final Solution
    for k = 2:n-1
        
        U1(k) = U1_old(k) + dU1dt(k)*dt;
        U2(k) = U2_old(k) + dU2dt(k)*dt;
        U3(k) = U3_old(k) + dU3dt(k)*dt;    
    end
    
    rho = U1./A;
    T = (gamma-1)*(U3./U1-0.5*gamma*(U2./U1).^2);
    v = U2./U1;
    
    
    % BCs
    
    % Inlet    
    U1(1) = A(1);
    U2(1) = 2*U2(2)-U2(3);
    U3(1) = U1(1)*(T(1)/(gamma-1)+0.5*gamma*(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);
    
    % Pressure
    p = rho.*T;
    
    % Mach Number 
    M = v./sqrt(T);
    
    % Mass Flow Rate
    m = U2;
    
    % Primitive Variables at throat 
    M_throat(i) = M(th);
    rho_throat(i) = rho(th);
    v_throat(i) = v(th);
    T_throat(i) = T(th);
    p_throat(i) = p(th);
    
    % Mass flow rate at different time steps
    figure(4)
    hold on
    if i ==1
        plot(x,m)
    end
    if i ==50
        plot(x,m)
    end
    if i ==100
        plot(x,m)
    end
    if i ==200
        plot(x,m)
    end
    if i ==400
        plot(x,m)
    end
    if i ==800
        plot(x,m)
    end
             
end

% Plotting 

% Mass flow rate at different time steps
xlabel('Nozzle Length')
ylabel('Mass Flow Rate')
legend('1 dt','50 dt','100 dt','200 dt','400 dt','800 dt')
title('Variation of Mass Flow Rate at different time steps')
% Steady State Primitive Variable Profiles

figure(1)
subplot(5,1,1)
plot(x,M,'b')
xlabel('Nozzle Length')
ylabel('Mach Number')
title('Variation of Steady State Primitive Variable Profiles')
subplot(5,1,2)
plot(x,rho,'k')
xlabel('Nozzle Length')
ylabel('Density')
subplot(5,1,3)
plot(x,v,'m')
xlabel('Nozzle Length')
ylabel('velocity')
subplot(5,1,4)
plot(x,T,'r')
xlabel('Nozzle Length')
ylabel('Temperature')
subplot(5,1,5)
plot(x,p,'y')
xlabel('Nozzle Length')
ylabel('Pressure')

% Transient State Primitive Variable Profiles
figure(2)
subplot(5,1,1)
plot(M_throat)
xlabel('Time Step')
ylabel('Mach Number')
title('Variation of Transient State Primitive Variable Profiles at Nozzle Throat')
subplot(5,1,2)
plot(rho_throat,'k')
xlabel('Time Step')
ylabel('Density')
subplot(5,1,3)
plot(v_throat,'r')
xlabel('Time Step')
ylabel('velocity')
subplot(5,1,4)
plot(T_throat,'g')
xlabel('Time Step')
ylabel('Temperature')
subplot(5,1,5)
plot(p_throat,'y')
xlabel('Time Step')
ylabel('Pressure')

% % Mass Flow rate
figure(3)
plot(x,m)
xlabel('Nozzle Length')
ylabel('Mass Flow Rate')
title('Mass Flow Rate along Nozzle length')
grid on

OUTPUT OF CONSERVATIVE FORM:

MATLAB CODE FOR COMPARISON OF CONSERVATIVE AND NON-CONSERVATIVE

FUNCTION CODE FOR NON-CONSERVATIVE

function [p,M,m] = non_conservative_form(x,dx,n,gamma,nt,CFL,rho,T,v,A);

% Time Loop
for i = 1:nt
     
    dt(i) = min((CFL*dx)./(v+sqrt(T)));
    
    rho_old = rho;
    T_old = T;
    v_old = v;
       
    % Predictive Method
    for j = 2:n-1
        
        dvdx = (v(j+1)-v(j))/dx;
        dlogA_dx = (log(A(j+1))-log(A(j)))/dx;
        drhodx  = (rho(j+1)-rho(j))/dx;
        dTdx = (T(j+1)-T(j))/dx;
        
        %Prdictive step
        
        % Continuity Equation
        drhodt_p(j) = -rho(j)*dvdx-rho(j)*v(j)*dlogA_dx-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)*dlogA_dx);
     
        %Predictive Solution update
        rho(j) = rho_old(j) + drhodt_p(j)*dt(i);
        v(j) = v_old(j) + dvdt_p(j)*dt(i);
        T(j) = T_old(j) + dTdt_p(j)*dt(i);
    end
    
    % Corrective Method
    
    for j = 2:n-1
        
        dvdx = (v(j)-v(j-1))/dx;
        dlogA_dx = (log(A(j))-log(A(j-1)))/dx;
        drhodx  = (rho(j)-rho(j-1))/dx;
        dTdx = (T(j)-T(j-1))/dx;
        
        %Prdictive step
        
        % Continuity Equation
        drhodt_c(j) = -rho(j)*dvdx-rho(j)*v(j)*dlogA_dx-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)*dlogA_dx);
     
        end
   
    % Average Values
    
    drhodt = 0.5*(drhodt_p+drhodt_c);
    dvdt = 0.5*(dvdt_p+dvdt_c);
    dTdt = 0.5*(dTdt_p+dTdt_c);
    
    % Final Solution
    for k = 2:n-1
        
        rho(k) = rho_old(k) + drhodt(k)*dt(i);
        v(k) = v_old(k) + dvdt(k)*dt(i);
        T(k) = T_old(k) + dTdt(k)*dt(i);    
    end
    
    % BCs
    % Inlet
    
    rho(1) = 1;
    T(1) = 1;
    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);
    
    % Pressure
    p = rho.*T;
    
    % Mach Number 
    M = v./sqrt(T);
    
    % Mass Flow Rate
    m = rho.*A.*v;
    
       
    
end
end

FUNCTION CODE FOR CONSERVATIVE

function [p,M,m] = conservative_form(x,dx,n,gamma,nt,CFL,A);

  for i = 1:n
    if x(i) <= 0.5 
    rho(i) = 1;
    T(i) = 1;
    end
    if 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);
    end
    if x(i) >=1.5
    rho(i) = 0.634-0.387*(x(i)-1.5);
    T(i) = 0.833-0.35*(x(i)-1.5);
    end
  end

 v = 0.579./(rho.*A);


U1 = rho.*A;
U2 = rho.*A.*v;
U3 = U1.*(T/(gamma-1)+0.5*gamma*v.^2);

% Time Loop
for i = 1:nt
     
    dt(i) = min((CFL*dx)./(v+sqrt(T)));
    
    U1_old = U1;
    U2_old = U2;
    U3_old = U3;
    
    F1=U2;
    F2=((U2.^2)./U1)+((gamma-1)/gamma)*(U3-(gamma/2)*(((U2).^2)./U1));
    F3=(((gamma)*(U2.*U3))./U1)-(((gamma*(gamma-1))/2)*((U2.^3)./U1.^2));
       
    % Predictive Method
    for j = 2:n-1
        
             
       dU1dt_p(j) = -(F1(j+1)-F1(j))/dx;
        
       J2 = (1/gamma)*rho(j)*T(j)*(A(j+1)-A(j))/dx;
       
       dU2dt_p(j) = -(F2(j+1)-F2(j))/dx+J2;
       
       dU3dt_p(j) = -(F3(j+1)-F3(j))/dx;
           
        %Predictive Solution update
        U1(j) = U1(j) + dU1dt_p(j)*dt(i);
        U2(j) = U2(j) + dU2dt_p(j)*dt(i);
        U3(j) = U3(j) + dU3dt_p(j)*dt(i);
        
    end
    
    %Update F & Pritimitive variables vector for J2
    rho = U1./A;
    T = (gamma-1)*(U3./U1-0.5*gamma*(U2./U1).^2);
    
    F1=U2;
    F2=((U2.^2)./U1)+((gamma-1)/gamma)*(U3-(gamma/2)*(((U2).^2)./U1));
    F3=(((gamma)*(U2.*U3))./U1)-(((gamma*(gamma-1))/2)*((U2.^3)./U1.^2));
    
% Corrective Method
    
    
    for j = 2:n-1        
             
             
       dU1dt_c(j) = -(F1(j)-F1(j-1))/dx;
        
       J2 = (1/gamma)*rho(j)*T(j)*(A(j)-A(j-1))/dx;
       
       dU2dt_c(j) = -(F2(j)-F2(j-1))/dx+J2;
       
       dU3dt_c(j) = -(F3(j)-F3(j-1))/dx;
        
            
        end
   
    % Average Values
    
    dU1dt = 0.5*( dU1dt_p+ dU1dt_c);
    dU2dt = 0.5*( dU2dt_p+ dU2dt_c);
    dU3dt = 0.5*( dU3dt_p+ dU3dt_c);
    
    % Final Solution
    for k = 2:n-1
        
        U1(k) = U1_old(k) + dU1dt(k)*dt(i);
        U2(k) = U2_old(k) + dU2dt(k)*dt(i);
        U3(k) = U3_old(k) + dU3dt(k)*dt(i);    
    end
    
    % BCs
    
    % Inlet    
    U1(1) = A(1);
    U2(1) = 2*U2(2)-U2(3);
    v(1)= U2(1)/U1(1);
    U3(1) = U1(1)*(T(1)/(gamma-1)+0.5*gamma*(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);
    
    rho = U1./A;
    T = (gamma-1)*(U3./U1- 0.5*gamma*(U2./U1).^2);
    v = U2./U1;
    
    % Pressure
    p = rho.*T;
    
    % Mach Number 
    M = v./sqrt(T);
    
    % Mass Flow Rate
    m = rho.*A.*v;
    
    
end
             
end

MAIN SCRIPT

clc 
clear all
close all

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

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

% Time Step
nt = 1400;
CFL = 0.5;

% Non Conservative Function call
[p_NC,M_NC,m_NC] = non_conservative_form(x,dx,n,gamma,nt,CFL,rho,T,v,A);

% Conservative Function call
[p_C,M_C,m_C] = conservative_form(x,dx,n,gamma,nt,CFL,A);

figure(1)
line([0 3],[0.579 0.579],'color','g')
hold on
plot(x,m_NC)
hold on 
plot(x,m_C)
xlabel('Nozzle Length')
ylabel('Mass Flow Rate')
legend('Analytical','Non Conservative','Conservative')
title('Mass Flow Rate along Nozzle length')
grid on

OUTPUT

Grid independence test :-

31 grid points

(n=16)

31 grid points

(n=16)

61 grid points

(n=31)

61 grid points

(n=31)

Flow Variables

at Throat

Analytical Solution

Conservative Form

Non-Conservative

Form

Conservative

Form

Non-Conservative

Form

Density 0.634 0.6524 0.6387 0.6390 0.6377
Temperature 0.833 0.8411 0.8365 0.8357 0.8354
Pressure 0.528 0.5487 0.5342 0.5340 0.5327
Mach Number

1.000

.9780 .9994 0.9978 0.9998

Above are the values for Density, Temperature, Pressure and Mach Number at the throat in analytical form, conservative and non conservative forms.

  • It is obsered that the non conservation form results are really closser to the analytical solution.
  • Also, on doubling the node points, the stedy state numerical results are tending toward the analytical solution.
  • Conservative numerical results are more closer to analytical solutions compared to non-conservative approach.

RESULTS/DISCUSSION

  • From the steeady state graphs of Distribution of primitive variables inside the nozzle ,
    • it is observed that as we move from the inlet to outlet the mach number goes on increaseing. At the inlet Mach number will be less than 1, at the throat the mach number will be equal to 1 and at the outlet the mach number will be greater than 1.
    • At the convergeing section the velocity will be less but velocity above mach number equal to 1 increases as the pipe cross section increases.
    • Energy is created from the heat and pressure of the fluid, as velocity increases the temperature and pressure decreases. As the pressure decreases the density also decreases.
  • From the transient state graphs of Distribution of primitive variables inside the nozzle
    • Initially the primitive values ocillates, after which the final,steady state value is approched.
    • These oscillations are associated with various unsteady compression and expansion waves which propagate through the nozzle during the transient process.
    • Later the time derivatives rapidly grow small over a finite number of timesteps.
  • From the graph of comparison between conservative and non conservative,
    • It is observed that conservation form of the equation is much more satisfactory than that obtained from the non conservation form.
    • Conservative numerical results are more closer to analytical solutions compared to non-conservative approach.

CONCLUSIONS

  • Conservative form is better compared to non-conservative.
  • For accuracy, non-conservative is better due to the direct calculation of primitive variables. While in conservation form primitive variables are calculated indirectly from flux variables. Increased number of calculations will lead to more errors. This is the case when we take primitive valrables. But if we are taking Flux variables conservation form is better.
  • The primitive variables converge to a numerical steady-state faster when using the non-conservation form.Because the dependent variables are primitive in non-conservation forms but not in conservation form, this might be the cause for that. Thus simulation time for the non-conservation form is shorter than the conservation form.
Your portfolio, next

Ship work like this
in your own course.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.