Simulation of 1D Super-sonic nozzle flow using MacCormack's Method
Objective
- To write code in MATLAB to solve and simulate the Quasi 1D Supersonic Nozzle Flow Equations using the MacCormack Method, by implementing both the conservative and non-conservative forms of the governing equations.
- To perform grid dependence test.
- To determine the minimum number of runs for which the simulation should be performed in order to attain convergence.
- To compare the normalized mass flow rate in the conservative and non-conservative forms.
- To determine which method is faster.
Subsonic - Supersonic Nozzle
In this problem, we consider the steady, isentropic (constant entropy) flow through a convergent-divergent nozzle as shown here. The inlet of the nozzle is considered as a reservoir of infinite area. Since area is infinite `(A rarr oo )` the velocity is very small `( V rarr 0)`. The flow expands isentropically to supersonic speeds at the nozzle exit. Thus, the flow is locally subsonic in the convergent section of the nozzle, sonic at the throat, and supersonic at the divergent section of the nozzle. The sonic flow (M=1) at the throat means that the local velocity at this section is equal to the local speed of sound. Thus, the given nozzle is called a subsonic-supersonic nozzle.
For the given nozzle, although the flow field is two dimensional in reality, we make an assumption that the flow properties vary only with x. Such flow is termed as a Quasi one dimensional flow (Quasi 1-D).
Problem Setup
In order to simulate the above problem numerically, the following steps are followed.
- The governing equations for the fluid flow are determined in terms of partial differential equations in order to perform a time marching process and obtain a time marching solution.
- The governing equations are discretized by applying a finite difference technique called the MacCormack's Method, and the corresponding expressions are solved.
- Other details such as time step calculations and boundary conditions setup are performed.
Governing Flow Equations
The governing equations for the given fluid flow problem are derived from the Navier-Stokes Equations. These governing equations are partial differential equations and are derived in both the conservative and non-conservative forms. The three main governing equations required for simulating the fluid flow are:
- The Continuity Equation
- The Momentum Equation
- The Energy Equation
These equations, in their conservative and non-conservative forms are presented below.
Governing Flow Equations - Non-conservative Form
1. Continuity Equation
`(delrho')/(delt') = -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')]`
where,
`T' = T/T_o`
`rho'=rho/rho_o`
`x'=x/L`
`V' = V/a_o`
`t' = t/(L/a_o)`
`A' = A/A^**`
The above equations have been take from "Computational Fluid Dynamics: Basics with Applications" by John D Anderson. The flow field variables have been non-dimensionalized accordingly.
The conservation forms of the governing equations are given as follows.
Governing Flow Equations - Conservation Form
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'((e')/(gamma-1)+gamma/2V'^2)A'])/(delt') + (del[rho'((e')/(gamma-1)+gamma/2V'^2)V'A'+p'A'V'])/(delx') = 0`
The above governing equations are in the conserved form and have been non-dimensionalized accordingly, as per the earlier equations.
To express these governing equations in a more generic form, we replace the flow field variables as below.
`U_1 = rho'A'`
`U_2 = rho'A'V'`
`U_3 = rho'(e/(gamma-1) + gamma/2V'^2)A'`
`F_1 = rho'A'V'`
`F_2 = rho'A'V'^2 + 1/gamma p'A'`
`F_3 = rho'((e')/(gamma-1)+gamma/2V'^2)V'A'+p'A'V'`
`J_2 = 1/gamma p'(delA')/(delx')`
Thus, by substituting the above equations in the governing equations, we get
`(delU_1)/(delt') = -(delF_1)/(delx')` - Continuity Equation
`(delU_2)/(delt') = -(delF_2)/(delx')+J_2` - Momentum Equation
`(delU_3)/(delt') = -(delF_3)/(delx')` - Energy Equation
The above governing flow equations are discretized by an appropriate finite difference technique (MacCormack's Method) to perform a time marching process and obtain a time-dependent solution to the flow parameters.
MacCormack's Method
The MacCormack's Method is an explicit finite difference technique that is second order accurate in both space and time. MacCormack's method is a widely used discretization scheme for obtaining a time marching solution of hyperbolic partial differential equation. The MacCormack's method is also called as the Predictor-Corrector Method, as it contains two steps towards evaluating the finite difference expressions, namelu the Predictor Step and the Corrector Step. The two steps of the MacCormack's method are described below.
Consider a one dimensional grid containing 'n' points. Let the flow variable to be evaluated be the flow velocity 'u'.
According to the MacCormack's Method, the value of the flow variable (u) at the next time step (`t+Deltat`) using the value at the current time step (t), is evaluated as,
`u_i^(t+Deltat) = u_i^t + ((delu)/(delt))_(avg) .Deltat`
where, `((delu)/(delt))_(avg)` is a representative of the mean value of `(delu)/(delt)` between the times t and `t+Deltat`
The average value of `((delu)/(delt))_(avg)` is given as,
`((delu)/(delt))_(avg) = 0.5(((delu)/(delt))_p+((delu)/(delt))_c)`
`((delu)/(delt))_p` = Derivative obtained from the Predictor Step
`((delu)/(delt))_c` = Derivative obtained from the Corrector Step
To describe the Predictor and Corrector Steps, consider a flow variable 'f'. Here, 'f' can denote anything (velocity, density, temperature etc.)
Predictor Step
To evaluate the Predictor step derivative, the spatial derivatives in the equations are replaced by their respective forward difference expressions as follows.
`((delf)/(delt))_p^t = (f_(i+1)^t - f_(i)^t)/(Deltax)`
All flow variables at time 't' are known i.e. the right hand side is known. Now a predicted value of the flow variable 'f' at the next time step `t+Deltat` is obtained using the Taylor's series expansion as follows.
`f_i^(t+Deltat) = f_i^t + ((delf)/(delt))_p^t* Deltat`
The value obtained above is a predicted value and is only first order accurate.
Corrector Step
In the corrector step, we obtain the predicted value of the time derivative at time `t+Deltat` by substituting the predicted value of the flow variable 'f' into the right side, replacing the spatial derivatives with backward difference expressions. (Updating the values and finding the derivative)
`((delf)/(delt))_c^(t+Deltat) = (f_(i)^(t+Deltat)-f_(i-1)^(t+Deltat))/(Deltax)`
The average value of the time derivative of the flow variable is the arithmetic mean of the time derivatives obtained from the predictor step and corrector step.
`((delf)/(delt))_(avg) = 0.5*(((delf)/(delt))_p+((delf)/(delt))_c)`
This allows us to obtain the "corrected' value of the flow variable at time `t+Deltat` using the equation given below.
`f_i^(t+Deltat) = f_i^t+((delf)/(delt))_(avg).Deltat`
The predictor-corrector sequence describe above yields the value of the flow variable 'f' at the grid point 'i' at the time step '`t+Deltat`'. This technique is used to calculate all the flow variables such as velocity, density, Temperature etc.
The MacCormack's Method as described above is a second order accurate method, as it uses a two step predictor-corrector sequencing utilizing both forward and differencing approaches to obtain the finite difference expressions. This method can be applied very easily as it does not require evaluating the second order derivatives of the flow variables directly. Hence, the governing equations in the conservative and non-conservative form are discretized in a one dimensional grid, using the MacCormack's Method, in order to peform a time marching apporach and obtain the flow parameters at each time step, and at all the grid points.
Time Step Calculation
The governing equations as described above are hyperbolic partial differential equations with respect to time. Hence, a stability constraint exists on this system, which is given as,
`Delta t = C((Deltax)/(a+V))`
where, C is the Courant Number
a is the local velocity of sound
V is the flow velocity
`Deltax` is the grid spacing
From the stability analysis criterial given by the CFL Condition, for a linear hyperbolic equation,
`C <=1` in order to obtain a stable numerical solution.
For the present case, the Courant Number C is assumed to be 0.5 to calculate the value of dt.
Boundary Conditions
- Inflow Boundary Conditions
The boundary conditions at the nozzle inlet is given as follows
`V_1 = 2V_2 -V_3`
- Outflow Boundary Conditions
The boundary conditions at the nozzle outlet are given as follows
`V_n = 2V_(n-1) - V_(n-2)`
`rho_n = 2rho_(n-1) - rho_(n-2)`
`T_(n) = 2T_(n-1) -T_(n-2)`
Nozzle Shape and Initial Conditions
The Nozzle Shape A = A(x) is specified and held fixed, independent of time. For the given case, a parabolic area distribution is assumed as follows.
`A = 1+2.2(x-1.5)^2` `0<=x<=3`
`x<1.5` - Convergent Portion
`x>1.5` - Divergent Portion
`x = 1.5`- Throat
Initial conditions at time t=0 are given as,
`rho = 1-0.3146x`
`T = 1-0.2314x`
`V = (0.1+1.09x)T^(0.5)`
Using the initial conditions and boundary conditions described above, the MATLAB Code is written to solve the governing equations using the MacCormack's Method and obtain the solutions for the flow variables at each time step.
MATLAB Code for Main Program
clear all
close all
clc
%Quasi 1D Supersonic Flow through a Nozzle using MacCormack Method
n=31;
%n1 = 61 - for grid dependency test
x=linspace(0,3,n);
dx = x(2)-x(1); %Grid spacing
a = 1+(2.2)*(x-1.5).^2; %Area
throat = (n+1)/2;
C=0.5; %Courant Number
%No. of timesteps
nt = 1400;
gamma=1.4;
m_ana=0.579*ones(1,n); %Analytical Solution of Mass Flow
%Non - Conservative Form Solution
tic;
[rho_nc,T_nc,v_nc,M_nc,m_nc,P_nc,rho_tnc,v_tnc,T_tnc,P_tnc,M_tnc,m_tnc]=noncon(n,x,dx,nt,C,gamma,a,throat);
nc_simtime=toc;
%Conservative Form Solution
tic;
[rho_c,T_c,v_c,m_c,M_c,P_c,rho_tc,v_tc,T_tc,P_tc,M_tc,m_tc]=conserv(n,x,dx,nt,C,gamma,a,throat);
c_simtime=toc;
%Plotting variation of flow variables across nozzle length - Non conservative form
figure(3)
plot(x,T_nc,'linewidth',2.5);
hold on
plot(x,T_c,'*');
grid on
xlabel('Nozzle Length');
ylabel('Temperature');
title('Variation of Temperature across Nozzle Length');
legend('Non-conservative','Conservative');
figure(4)
plot(x,M_nc,'linewidth',2.5);
hold on
grid on
plot(x,M_c,'*');
xlabel('Nozzle Length');
ylabel('Mach Number');
title('Variation of Mach Number across Nozzle Length');
legend('Non-conservative','Conservative');
figure(5)
plot(x,v_nc,'linewidth',2.5);
hold on
grid on
plot(x,v_c,'*');
xlabel('Nozzle Length');
ylabel('Velocity');
title('Variation of Velocity across Nozzle Length');
legend('Non-conservative','Conservative');
figure(6)
plot(x,P_nc,'linewidth',2.5);
hold on
grid on
plot(x,P_c,'*');
xlabel('Nozzle Length');
ylabel('Pressure');
title('Variation of Pressure across Nozzle Length');
legend('Non-conservative','Conservative');
figure(7)
plot(x,rho_nc,'linewidth',2.5);
hold on
grid on
plot(x,rho_c,'*');
xlabel('Nozzle Length');
ylabel('Density');
legend('Non-conservative','Conservative');
title('Variation of Density across Nozzle Length');
%Plotting the variation of Flow variables across the throat - Non conservative form
ntgrid=linspace(1,1400,1400);
figure(8)
subplot(3,2,1)
plot(ntgrid,T_tnc);
grid on
xlabel('Nozzle Length');
ylabel('Temperature');
title('Temperature variation across throat');
subplot(3,2,2)
plot(ntgrid,M_tnc);
grid on
xlabel('Nozzle Length');
ylabel('Mach Number');
title('Mach Number variation across throat');
subplot(3,2,3)
plot(ntgrid,v_tnc);
grid on
xlabel('Nozzle Length');
ylabel('Velocity');
title('Velocity variation across throat');
subplot(3,2,4)
plot(ntgrid,P_tnc);
grid on
xlabel('Nozzle Length');
ylabel('Pressure');
title('Pressure variation across throat');
subplot(3,2,5)
plot(ntgrid,rho_tnc);
grid on
xlabel('Nozzle Length');
ylabel('Density');
title('Density variation across throat');
subplot(3,2,6)
plot(ntgrid,m_tnc);
grid on
xlabel('Nozzle Length');
ylabel('Mass Flow');
title('Mass flow variation across throat');
%Plotting the variation of Flow variables across the throat - Conservative form
ntgrid=linspace(1,1400,1400);
figure(9)
title('Mass flow variation across throat (Conservative)');
subplot(3,2,1)
plot(ntgrid,T_tc);
grid on
xlabel('Nozzle Length');
ylabel('Temperature');
title('Temperature variation across throat');
subplot(3,2,2)
plot(ntgrid,M_tc);
grid on
xlabel('Nozzle Length');
ylabel('Mach Number');
title('Mach Number variation across throat');
subplot(3,2,3)
plot(ntgrid,v_tc);
grid on
xlabel('Nozzle Length');
ylabel('Velocity');
title('Velocity variation across throat');
subplot(3,2,4)
plot(ntgrid,P_tc);
grid on
xlabel('Nozzle Length');
ylabel('Pressure');
title('Pressure variation across throat');
subplot(3,2,5)
plot(ntgrid,rho_tc);
grid on
xlabel('Nozzle Length');
ylabel('Density');
title('Density variation across throat');
subplot(3,2,6)
plot(ntgrid,m_tc);
grid on
xlabel('Nozzle Length');
ylabel('Mass Flow');
%Comparison of conservative, non-conservative and analytical solutions of mass flow
figure(10)
plot(x,m_c,'linewidth',2.5);
hold on
grid on
plot(x,m_nc,'linewidth',2.5);
hold on
plot(x,m_ana,'*');
xlabel('Nozzle Length');
ylabel('Mass Flow');
legend('Conservative Form Solution','Non-conservative Form Solution','Analytical Solution');
title('Mass Flow rate comparison of conservative, non-conservative, analytical vs Nozzle Length');
MATLAB Code for Function to solve Non-conservative form of the governing equations
function[rho_nc,T_nc,v_nc,M_nc,m_nc,P_nc,rho_tnc,v_tnc,T_tnc,P_tnc,M_tnc,m_tnc]=noncon(n,x,dx,nt,C,gamma,a,throat)
%Calculate initial profiles
%Profiles are assumed for velocity, density, area and pressure
rho_nc = 1-(0.3146*x); %Density
T_nc = 1-(0.2314*x); %Temperature
v_nc = (0.1+1.90*x).*T_nc.^0.5; %Velocity
m_ana=0.579*ones(1,n); %Analytical Solution for Mass Flow
m_initial=rho_nc.*a.*v_nc;
M_nc=v_nc./(T_nc.^0.5);
C=0.5; %Courant Number
for i=1:n
ts(i) = C*dx/((T_nc(i)^0.5)+v_nc(i));
end
dt = min(ts);
%outer time loop
for k=1:nt
rho_old=rho_nc;
v_old=v_nc;
T_old=T_nc;
%predictor method (forward differencing)
for j=2:(n-1)
dvdx = (v_nc(j+1)-v_nc(j))/dx;
dlogadx = (log(a(j+1))-log(a(j)))/dx;
drhodx = (rho_nc(j+1)-rho_nc(j))/dx;
dTdx = (T_nc(j+1)-T_nc(j))/dx;
%continuity equation
drhodt_p(j) = (-rho_nc(j)*dvdx)-(rho_nc(j)*v_nc(j)*dlogadx)-(v_nc(j)*drhodx);
%Momentum Equation
dvdt_p(j) = -v_nc(j)*dvdx-(1/gamma)*(dTdx + (T_nc(j)/rho_nc(j))*drhodx);
%Energy Equation
dTdt_p(j) = -v_nc(j)*dTdx-((gamma-1)*T_nc(j)*(dvdx+(v_nc(j)*dlogadx)));
%Solution Update
rho_nc(j) = rho_nc(j) + drhodt_p(j)*dt;
v_nc(j) = v_nc(j) + dvdt_p(j)*dt;
T_nc(j) = T_nc(j) + dTdt_p(j)*dt;
end
%Corrector Method (backward differencing)
for j=2:(n-1)
dvdx = (v_nc(j)-v_nc(j-1))/dx;
dlogadx = (log(a(j))-log(a(j-1)))/dx;
drhodx = (rho_nc(j)-rho_nc(j-1))/dx;
dTdx = (T_nc(j)-T_nc(j-1))/dx;
%continuity equation
drhodt_c(j) = (-rho_nc(j)*dvdx)-(rho_nc(j)*v_nc(j)*dlogadx)-(v_nc(j)*drhodx);
%Momentum Equation
dvdt_c(j) = -v_nc(j)*dvdx-(1/gamma)*(dTdx + (T_nc(j)/rho_nc(j))*drhodx);
%Energy Equation
dTdt_c(j) = -v_nc(j)*dTdx-((gamma-1)*T_nc(j)*(dvdx+(v_nc(j)*dlogadx)));
%Solution update need not be done for Corrector Method
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 i=2:n-1
v_nc(i) = v_old(i)+dvdt(i)*dt;
rho_nc(i) = rho_old(i) +drhodt(i)*dt;
T_nc(i) = T_old(i)+dTdt(i)*dt;
end
%Apply Boundary conditions
%Inlet BC
v_nc(1) = 2*v_nc(2) - v_nc(3);
%Outlet BC
v_nc(n) = 2*v_nc(n-1) - v_nc(n-2);
rho_nc(n) = 2*rho_nc(n-1) - rho_nc(n-2);
T_nc(n) = 2*T_nc(n-1) - T_nc(n-2);
P_nc=rho_nc.*T_nc;
M_nc=v_nc./T_nc.^0.5;
m_nc=rho_nc.*a.*v_nc;
%Flow variables at Throat
rho_tnc(k) = rho_nc(throat);
v_tnc(k) = v_nc(throat);
T_tnc(k)=T_nc(throat);
M_tnc(k)=M_nc(throat);
m_tnc(k)=m_nc(throat);
P_tnc(k)=P_nc(throat);
%Variation of mass at different time steps
figure(1)
if(k==10)
plot(x,m_nc,'r');
hold on
elseif(k==50)
plot(x,m_nc,'b');
hold on
grid on
elseif(k==100)
plot(x,m_nc,'g');
hold on
elseif(k==250)
plot(x,m_nc,'y');
hold on
elseif(k==500)
plot(x,m_nc,'k');
hold on
elseif(k==1000)
plot(x,m_nc,'m');
hold on
elseif(k==1400)
plot(x,m_nc,'c');
hold on
plot(x,m_ana,'--','linewidth',2.5);
hold on
grid on
title('Variation of Mass flow at different time steps - Non conservative form');
xlabel('Length of Nozzle');
ylabel('Mass Flow');
legend('10th Time Step','50th Time Step','100th Time Step','250th Time Step','500th Time Step','1000th Time Step','1400th Time Step','Analytical Solution');
end
end
end
MATLAB Code for Function to solve Conservative form of the governing equations
function[rho_c,T_c,v_c,m_c,M_c,P_c,rho_tc,v_tc,T_tc,P_tc,M_tc,m_tc] = conserv(n,x,dx,nt,C,gamma,a,throat)
%Calculate initial profiles
%Profiles are assumed for velocity, density, area and pressure
rho = zeros(1,n);
T = zeros(1,n);
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))
rho(i)=0.634-(0.3879*(x(i)-1.5));
T(i)=0.833-(0.3507*(x(i)-1.5));
end
end
%Pressure
P = rho.*T;
%Velocity
%rho.V.a = 0.59;
%V=0.59/rho.a;
v = 0.59./(rho.*a);
%Mass Flow
m = rho.*a.*v;
m_ana=0.579*ones(1,n); %Analytical Solution of Mass Flow
%Mach Number
M = v./(T.^0.5);
U1 = rho.*a;
U2 = U1.*v;
U3 = U1.*((T./(gamma-1))+((gamma/2).*(v.^2)));
for i=1:n
ts(i) = (C*dx)/(T(i)^0.5+v(i));
end
dt = min(ts);
%Predictor-Collector calculation
%Outer time loop
for k=1:nt
U1old=U1;
U2old=U2;
U3old=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));
%Predictor Method - Forward differencing
for j=2:(n-1)
df1dx_p = (F1(j+1)-F1(j))/dx;
df2dx_p = (F2(j+1)-F2(j))/dx;
df3dx_p = (F3(j+1)-F3(j))/dx;
dadx_p = (a(j+1)-a(j))/dx;
J2(j) = (1/gamma)*((rho(j).*T(j))*dadx_p);
%Continuity Equation
du1dt_p(j) = -df1dx_p;
%Momentum Equation
du2dt_p(j) = -df2dx_p + J2(j);
%Energy Equation
du3dt_p(j) = -df3dx_p;
%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
%Extracting Flow Variables
rho = U1./a;
v = U2./U1;
T = (gamma-1).*((U3./U1)-((gamma/2).*v.^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));
%Corrector Step - Backward Differencing
for j=2:(n-1)
df1dx_c = (F1(j)-F1(j-1))/dx;
df2dx_c = (F2(j)-F2(j-1))/dx;
df3dx_c = (F3(j)-F3(j-1))/dx;
dadx_c = (a(j)-a(j-1))/dx;
J2(j) = (1/gamma)*((rho(j).*T(j))*dadx_c);
%Continuity Equation
du1dt_c(j) = -df1dx_c;
%Momentum Equation
du2dt_c(j) = -df2dx_c + J2(j);
%Energy Equation
du3dt_c(j) = -df3dx_c;
end
%Solution Update not needed for Corrector Step
%Compute Average Time 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 i=2:(n-1)
U1(i) = U1old(i)+(du1dt(i)*dt);
U2(i) = U2old(i)+(du2dt(i)*dt);
U3(i) = U3old(i)+(du3dt(i)*dt);
end
%Boundary Conditions
%Inlet Boundary Conditions
U2(1) = 2*U2(2)-U2(3);
%Outlet Boundary Conditions
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)-(0.5*gamma*(v.^2)));
P = rho.*T;
M=v./(T.^0.5);
m=rho.*a.*v;
rho_c = rho;
v_c = v;
T_c = T;
P_c = P;
M_c = M;
m_c = m;
%Flow variables at Throat
rho_tc(k) = rho(throat);
v_tc(k) = v(throat);
T_tc(k) = T(throat);
M_tc(k) = M(throat);
m_tc(k) = m(throat);
P_tc(k) = P(throat);
%Variation of mass at different time steps
figure(2)
if(k==10)
plot(x,m_c,'r');
hold on
elseif(k==50)
plot(x,m_c,'b');
hold on
grid on
elseif(k==100)
plot(x,m_c,'g');
hold on
elseif(k==250)
plot(x,m_c,'y');
hold on
elseif(k==500)
plot(x,m_c,'k');
hold on
elseif(k==1000)
plot(x,m_c,'m');
hold on
elseif(k==1400)
plot(x,m_c,'c');
hold on
plot(x,m_ana,'--','linewidth',2.5);
hold on
grid on
title('Variation of Mass flow at different time steps - Conservative form');
xlabel('Length of Nozzle');
ylabel('Mass Flow');
legend('10th Time Step','50th Time Step','100th Time Step','250th Time Step','500th Time Step','1000th Time Step','1400th Time Step','Analytical Solution');
end
end
end
Plotting the Results
- Variation of Flow parameters over the Nozzle Length (Conservative vs. Non-conservative)
1. Variation of Temperature across Nozzle Length
2. Variation of Pressure across Nozzle Length
3. Variation of Density across Nozzle Length
4. Variation of Velocity across Nozzle Length
5. Variation of Mach Number across Nozzle Length
6. Variation of Mass Flow across Nozzle Length at different time steps
- Non-conservative Form
- Conservative Form
- Comparison of Steady State Mass distribution across Nozzle Length (Conservative vs. Non-conservative vs. Analytical Solutions)
- Simulation time for conservative and non-conservative forms
- Variation of Flow parameters across the Nozzle Throat at different time steps
Non-conservative Form
Conservative Form
Inference
- From the plots comparing the variation of flow paramters across the nozzle length, we can see that the conservative and non-conservative plots for Density, Velocity, Temperature, Mach Number, and Pressure are almost identical. This implies that the conservative and non-conservative forms of the governing equations are mathematically equivalent, and there are no differences in the results for this case.
- From the plots, it can also be seen that the calculated profiles for the flow parameters across Nozzle Length is in agreement with the expected flow profiles (i.e. Temperature Drop, Pressure Drop, Density Drop, Velocity increase, and Pressure increase)
- On comparing the conservative and non-conservative plots for the variation of mass flow across nozzle length at different time steps, it can be seen that, the non-conservative form plot has more oscillations during the initial time steps (50 and 100), when compared to the conservative form. These oscillations eventually subside at subsequent time steps leading to a converged solution. In addition to this, there are a lot of oscillations at the nozzle inlet and outlet for the non-conservative form.
- From the conservative and non conservative mass distribution plots, it can be seen that the non-conservative form converges close to the 1000th time step (approx. 900) whereas, the non-conservative form converges at around 700 time steps. The computational time for solving the conservative form is 1.97 secs, whereas the time for solving the non-conservative form is 2.02 secs. Thus, from these results, it can be concluded that the conservative form is a faster method than the non-conservative form.
- From the plot comparing the steady state mass distribution of conservative, non-conservative and the analytical solutions, it is evident that the conservative form solution is more closer to the analytical form than the non-conservative form. The non-conservative form also has a number of oscillations, particularly at the nozzle inlet and outlet, whereas the conservative form remains a lot more stable throughout the nozzle length. This is because, in the conservative form, the mass flow itself is one of the dependent variables in the governing equations, unlike the non-conservative form where the primitive flow variables (velocity, density) are the dependent variables, which are solved and then used to calculate the mass flow. In essence, the conservative form does a better job of conserving or preserving the mass flow than the non-conservative form.
- From the plot comparing the variation of the flow parameters across the nozzle throat, it can be seen for both the conservative and non-conservative forms that the flow variables oscillate during the initial time steps and then become stable at subsequent time steps.
- From the non-conservative plots of the flow variables, it can be seen that the flow variables across the throat attain steady state convergence at around 900 time steps. From the conservative plots, it can be seen that the flow variables across the throat attain steady state convergence at around 750-800 time steps. This implies, the convergence rate of conservative form is a lot more faster than the non-conservative form.
- From all the above plots,
Grid Dependency Test
The matter of grid dependency is a serious consideration in CFD. If the numerical solution obtained changes with a change in number of grid points, then the solution is grid dependent. In this case, the grid points should be increased until the solution no longer changes with the number of grid points. This leads to a grid independent solution, which is desirable.
In order to test the grid dependency of the numerical solution obtained, the no. of grid points are increased from 31 to 61, and the solution procedure is performed again utilizing both the conservative and non-conservative forms of the equations. The steady state results are compared in order to qualitatively evaluate the grid dependency of the numerical solution obtained.
Consider the steady state numerical solution of the flow variables at the throat, from both the conservative and non-conservative forms. Tabulating these values and comparing them with the analytical solution values, we get the following table.
Inference
- The above table represents the values of the flow variables from the conservative and non-conservative form solutions of the governing equations.
- From the above table, we can clearly see that the values for the flow parameters obtained from the non-conservative form of the governing equation, are closer to the analytical results than the conservative form of the equations, at 31 grid points.
- As the no. of grid points is increased to 61, there is a marginal change in the flow parameter values. The flow values from the conservative forms get more closer to the analytical solution, but they are not as close as the non-conservative form solutions, with half the grid points (31 points)
- Thus, it is clear that irrespective of the no. of grid points, the non-conservative form solution produces more accurate values to the flow parameters.
- In essence, the conservative form yields a better mass flow distribution, whereas the non-conservative form yields more accurate flow values that are closer to the analytical solution.