Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
Objective
To solve Quasi 1D Converging-Diverging Nozzle numerically for both conservative & non-conservative forms of governing equations by implementing MacCormack's Method using MATLAB.
Problem Description
We consider a steady, isentropic flow through the converging-diverging nozzle. The inlet of the nozzle is connected with the reservoir. The cross-sectional area of the reservoir is large (tending to infinity), and hence the velocity inside the reservoir is very small(tending to 0). The flow is locally subsonic in the converging section of the nozzle, sonic at the nozzle throat, and expands isentropically to the supersonic speed at the nozzle exit.
In this project, non-dimensional governing equations for the Quasi 1D Converging-Diverging Nozzle in the conservative and non-conservative form are discretized & solved with MacCormack's Method using MATLAB.
The following numerical results obtained from both forms of equations are compared:-
1. Non-Dimensional mass flow rate across the nozzle for different time steps.
2. Transient variation of flow variables at the nozzle's throat location.
3. Steady-state flow field across with the nozzle.
4. Comparison of steady-state mass flow rate along the nozzle x-direction for Conservative & Non-Conservative form.
5. Comparison of steady-state flow variables at the throat location between Conservative, Non-Conservative and Analytical results.
6. Grid sensitivity study.
Non-Dimensional Non-Conservative Form of governing equation
Continuity Equation
`(partialrho')/(partialt') = - rho'.(partialV')/(partialx')-rho'.V.'(partiallnA')/(partialx')-V'.(partialrho')/(partialx')`
Momemtum Equation
`(partialV')/(partialt') = - V'.(partialV')/(partialx')-1/gamma.((partialT')/(partialx')+(T')/(rho').(partialrho')/(partialx'))`
Energy Equation
`(partialT')/(partialt') = - V'.(partialT')/(partialx')-(gamma-1).T'.((partialV')/(partialx')+V'.(partiallnA')/(partialx'))`
Non-Dimensional Conservative Form of governing equation
Continuity Equation
`(del(rho'A'))/(delt') + (del(rho'A'V'))/(delx') = 0`
Momemtum Equation
`(del(rho'A'V'))/(delt') + (del[rho'A'(V')^2 +(1/gamma)P'A'])/(delx') = (1/gamma)P'((delA')/(delx'))`
Energy Equation
`(del[rho'((e')/(gamma-1) + (gamma/2)(V')^2)A'])/(delt') + (del[rho'((e')/(gamma-1) + (gamma/2)(V')^2)V'A' + P'A'V'])/(delx') = 0`
Gas constant for air, `gamma = 1.4`
MacCormack's Method
MacCormack's Method consists of two-step, First the Predictor step & then the Corrector step.
Let us consider a variable `f(x)` which can be any flow variable such as velocity, pressure, density, temperature, etc.
1. To start with MacCormack Method, First, we need the initial profile of the flow variables.
2. Predictor Step - Here we will calculate the time derivative using the initial profile of variables & the forward difference scheme.
`((delf)/(dt))_p = (f(i+1)_("Initial") - f(i)_("Initial"))/(Deltax)`
3. Solution update for the predictor step
`f_p = f_(Initial) + ((delf)/(delt))_p * Deltat`
`f_p` = value of flow variable at the predicted step
`f_("Initial")`= Intial value of the flow variable
`delt` = Time Step size.
4. Corrector Step - Here we will calculate the time derivative using the predicted profile of variables & the backward difference scheme.
`((delf)/(dt))_c = (f(i)_p - f(i-1)_p)/(Deltax)`
5. Average time derivative calculation.
`((delf)/(dt))_(avg) = 0.5(((delf)/(dt))_p +((delf)/(dt))_c )`
6. Final Solution update for the next time step.
`(f_i)^(n+1) = (f_(Initial))^n + ((delf)/(delt))_(avg) * Deltat`
`(f_i)^(n+1)` = value of solution variable at the next time step
Time Step Size Calculation
`delt = CFL((delx)/(a + v))`
`delt =` time step size
`delx`= minimum mesh size
a = local speed of sound inside the nozzle
v = local flow speed inside the nozzle.
`a = sqrt(T')`
`T'`= Local flow temperature inside the nozzle
`delt` will be calculated for each and every grid point, from i = 1 to i = N based on the local flow and sound speed and minimum value of time step size will be chosen for the next time marching step.
`delt = "minimum"(delt_1,delt_2,delt_3,delt_4....... delt_n)`
Boundary Condition
Unsteady, Inviscid flow is governed by hyperbolic equations, and therefore for 1D unsteady flow there exist two real characteristic lines through any point in xt plane. Physically these two characteristics represent infinitely week Mach waves that are propagating upstream and downstream. Both Mach waves are traveling at the speed of sound 'a'. Let's have a look at the image below where we have a converging-diverging nozzle with xt plane drawn below it. This image will help us to define inlet and outlet boundary conditions.
The method of characteristics tells us that at a boundary where one characteristic propagates into the domain, then the value of one variable must be specified at that boundary and if one characteristics line propagates out of the domain, then the value of another flow variable must be extrapolated from the inside domain during the time marching solution.
Inflow BC - If we look at point 1(Inlet of the nozzle), the local flow velocity is subsonic (v < a). The left running Mach wave, which is traveling toward the left(relative to fluid movement) at the speed of sound easily works its way upstream against the low subsonic flow, which is moving from left to right. Hence the left running characteristics are shown as `a_1 - v_1` (relative to the fixed nozzle). So the left running characteristics propagating out of the domain. At point 1, the right running characteristics is entering into the domain with speed `a_1 + v_1` because fluid is traveling towards the right and right running Mach(characteristics) wave is also traveling towards the right at speed of sound relative to the fluid. Also, note at point 1, Streamline is traveling from left to right, Streamline plays the same role as a characteristic line. Hence at point 1 Streamline is entering into the domain.
Hence at the inlet(point 1), We must the specific value of two dependent flow variables as one characteristic line and streamline entering into the domain and one flow variable must be extrapolated from the inside domain as one characteristic line going out from the inlet.
`rho_1' = 1` [Fixed value at the inlet]
`T_1' = 1` [Fixed value at the inlet]
`v_1=2*v_2 - v_3` [Extraolated from the inside domain]
velocity is chosen to be float at inlet because on a physical basis we know the mass flow through the nozzle must be allowed to adjust to the proper steady-state and allowing `V_1'` to adjust based on the flow field growth in time is most sensible.
Since point 1 is viewed as the reservoir, we stipulate density and temperature value at point 1 to be their respective stagnation value `rho_o` & `T_o`respectively.
Outflow BC - Now if we look at point N, Left running characteristics propagates left at the speed of sound 'a'. However the flow at the outlet is supersonic (v > a), the left running characteristics is carried downstream by the flow at a speed `v_n -a_n`, the right running characteristics propagates right at the speed of sound 'a' hence it is swept at a speed of `v_n + a_n`. So at the outlet both left and right running characteristics propagate out of the domain, so does the streamline at point N.
Hence at the outlet(point N), We must extrapolate all flow variables from the inside domain as both characteristic lines & streamline are going out from the outlet.
`v_n = 2v_(n-1) - v_(n-2)` [Extraolated from the inside domain]
`rho_n = 2rho_(n-1) - rho_(n-2)` [Extraolated from the inside domain]
`T_n = 2T_(n-1) -T_(n-2)` [Extraolated from the inside domain]
MATLAB code for plotting various numerical solution obtained by solving Conservative & Non-Conservative form of equations
clc
clear all
close all
%Rajat Walia
%SOLVING QUASI 1D SUPERSONIC NOZZLE FLOW
%All flow field variables are calculated and plotted in their non-dimensional form
%Inputs
n = 31; %mesh points
x = linspace(0,3,n); %1d mesh
dx = x(2) - x(1); %grid spacing
gamma = 1.4; %gas constant
nt = 800; %number of time step
CFL = 0.5; %CFL No.
% NON - CONSERVATIVE FORM
[v,pressure,T,Mach_no,mass_flow,rho,rho_throat,P_throat,v_throat,T_throat,Mach_no_throat, mass_flow_throat] = non_conservative_form(n,x,dx,gamma,CFL,nt);
%Potting time wise variation of non-dimensional primitive variables at the nozzle throat
hold on
figure(4)
%density
subplot(4,1,1)
plot(linspace(1,nt,nt),rho_throat,'color','m', 'LineWidth', 1.5)
ylabel('Density')
axis([0 1000 0 1.4]);
grid minor;
title("Non Dimensional Flow Variable variation with Time at Nozzle's Throat [Non Conservative Form]")
%temperature
subplot(4,1,2)
plot(linspace(1,nt,nt),T_throat,'color','c', 'LineWidth', 1.5)
ylabel('Temperature')
axis([0 1000 0.5 1]);
grid minor;
%pressure
subplot(4,1,3)
plot(linspace(1,nt,nt),P_throat,'color','g', 'LineWidth', 1.5')
ylabel('Pressure')
axis([0 1000 0.3 1.2]);
grid minor;
%Mach number
subplot(4,1,4)
plot(linspace(1,nt,nt),Mach_no_throat,'color','y', 'LineWidth', 1.5)
xlabel('Time Step')
ylabel('Mach No.')
axis([0 1000 0.5 1.5]);
grid minor;
%Potting steady-state field of non-dimensional primitive variables across nozzle
hold on
figure(5)
%density
subplot(4,1,1)
plot(x,rho,'color','m', 'LineWidth', 1.5)
ylabel('Density')
axis([0 3 0 1]);
grid minor;
title("Steady-State Flow Field vs Nozzle X-Direction [Non Conservative Form]")
%temperature
subplot(4,1,2)
plot(x,T,'color','c', 'LineWidth', 1.5)
ylabel('Temperature')
axis([0 3 0 1]);
grid minor;
%pressure
subplot(4,1,3)
plot(x,pressure,'color','g', 'LineWidth', 1.5')
ylabel('Pressure')
axis([0 3 0 1]);
grid minor;
%Mach number
subplot(4,1,4)
plot(x,Mach_no,'color','y', 'LineWidth', 1.5)
xlabel('Nozzle X-Direction')
ylabel('Mach No.')
axis([0 3 0 4]);
grid minor;
mass_flow_rate_NC = mass_flow;
T_NC = T_throat;
rho_NC = rho_throat;
P_NC = P_throat;
Mach_NC = Mach_no_throat;
%CONSERVATIVE - FORM
[v,pressure,T,Mach_no,mass_flow,rho,rho_throat,P_throat,v_throat,T_throat,Mach_no_throat, mass_flow_throat] = conservative_form(n,x,dx,gamma,CFL,nt);
%Potting time wise variation of non-dimensional primitive variables at the nozzle throat
hold on
figure(7)
%density
subplot(4,1,1)
plot(linspace(1,nt,nt),rho_throat,'color','m', 'LineWidth', 1.5)
ylabel('Density')
axis([0 1000 0 1.4]);
grid minor;
title("Non Dimensional Flow Variable variation with Time at Nozzle's Throat [Conservative Form]")
%temperature
subplot(4,1,2)
plot(linspace(1,nt,nt),T_throat,'color','c', 'LineWidth', 1.5)
ylabel('Temperature')
axis([0 1000 0.5 1]);
grid minor;
%pressure
subplot(4,1,3)
plot(linspace(1,nt,nt),P_throat,'color','g', 'LineWidth', 1.5')
ylabel('Pressure')
axis([0 1000 0.3 1.2]);
grid minor;
%Mach number
subplot(4,1,4)
plot(linspace(1,nt,nt),Mach_no_throat,'color','y', 'LineWidth', 1.5)
xlabel('Time Step')
ylabel('Mach No.')
axis([0 1000 0.5 1.5]);
grid minor;
%Potting steady-state field of non-dimensional primitive variables across the nozzle
hold on
figure(8)
%density
subplot(4,1,1)
plot(x,rho,'color','m', 'LineWidth', 1.5)
ylabel('Density')
axis([0 3 0 1]);
grid minor;
title("Steady-State Flow Field vs Nozzle X-Direction [Conservative Form]")
%temperature
subplot(4,1,2)
plot(x,T,'color','c', 'LineWidth', 1.5)
ylabel('Temperature')
axis([0 3 0 1]);
grid minor;
%pressure
subplot(4,1,3)
plot(x,pressure,'color','g', 'LineWidth', 1.5')
ylabel('Pressure')
axis([0 3 0 1]);
grid minor;
%Mach number
subplot(4,1,4)
plot(x,Mach_no,'color','y', 'LineWidth', 1.5)
xlabel('Nozzle X-Direction')
ylabel('Mach No.')
axis([0 3 0 4]);
grid minor;
mass_flow_rate_C = mass_flow;
T_C = T_throat;
rho_C = rho_throat;
P_C = P_throat;
Mach_C = Mach_no_throat;
%Comparison of normalized steady-state mass-flow rate across the nozzle
%between Conservative & Non-Conservative Form
%analytical non dimensional mass flow rate throught the nozzle = 0.590
mass_flow_rate_analytical = 0.590*ones(1,n);
figure(9)
hold on
plot(x, mass_flow_rate_NC, 'r', 'LineWidth', 1.5);
hold on
plot(x, mass_flow_rate_C, 'k', 'LineWidth', 1.5);
hold on
plot(x, mass_flow_rate_analytical, 'g', 'LineWidth', 1.5, 'LineStyle','--');
grid minor;
legend("Non Conservative Form", "Conservative Form", "Analytical");
xlabel("Nozzle X-Direction");
ylabel("Non Dimensional Mass Flow Rate");
title("Comparison of steady Normalized Mass Flow Rate for Conservative & Non-Conservative Form")Code for Non - Conservative Form of ID Quasi Super-Sonic Nozzle Flow
function [v,pressure,T,Mach_no,mass_flow,rho,rho_throat,P_throat,v_throat,T_throat,Mach_no_throat, mass_flow_throat] = non_conservative_form(n,x,dx,gamma,CFL,nt)
%initial profile
%all flow variables are in their non-dimensional form
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; %converging diverging area
mass_flow_initial = rho.*A.*v; %mass flow rate
throat = find(A == 1); %find throat location
total_sim_time = 0; %initial flow time
%solving the discretized equations (Non-Conservative Form)
%Maccormack Method
for k = 1:nt; %time marching
% Storing initial data
T_old = T;
v_old = v;
rho_old = rho;
% Calculating the time step
dt = min((CFL.*dx)./(sqrt(T) + v));
% Predictor method
for j = 2:n-1;
dvdx = (v(j+1) - v(j))/dx;
dlnAdx = (log(A(j+1)) - log(A(j)))/dx;
drhodx = (rho(j+1) - rho(j))/dx;
dTdx = (T(j+1) - T(j))/dx;
%continutiy equation
drhodt_p(j) = -rho(j)*dvdx - rho(j)*v(j)*dlnAdx - v(j)*drhodx;
%momemtum 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)*dlnAdx);
%solution update
v(j) = v(j) + dvdt_p(j)*dt;
rho(j) = rho(j) + drhodt_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;
dlnAdx = (log(A(j)) - log(A(j-1)))/dx;
drhodx = (rho(j) - rho(j-1))/dx;
dTdx = (T(j) - T(j-1))/dx;
%continutiy equation
drhodt_c(j) = -rho(j)*dvdx - rho(j)*v(j)*dlnAdx - v(j)*drhodx;
%momemtum 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)*dlnAdx);
end
% compute the average time derivative
dvdt = 0.5*(dvdt_p + dvdt_c);
drhodt = 0.5*(drhodt_p + drhodt_c);
dTdt = 0.5*(dTdt_p + dTdt_c);
% Solution update
for i = 2:n-1
v(i) = v_old(i) + dvdt(i)*dt;
rho(i) = rho_old(i) + drhodt(i)*dt;
T(i) = T_old(i) + dTdt(i)*dt;
end
%boundary condition
%inlet
rho(1) = 1; %fixed value
T(1) = 1; %fixed value
v(1) = 2*v(2) - v(3); %extrapolated
%outlet
%extrapolated from the inside domain
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 the non-dimensional Mass Flow Rate, Pressure & Mach Number
mass_flow = rho.*A.*v;
pressure = rho.*T;
Mach_no = v./sqrt(T);
%calculating non-dimensional variables at the nozzle throat
rho_throat(k) = rho(throat);
v_throat(k) = v(throat);
T_throat(k) = T(throat);
P_throat(k) = pressure(throat);
mass_flow_throat(k) = mass_flow(throat);
Mach_no_throat(k) = Mach_no(throat);
total_sim_time = total_sim_time + dt; %updating total flow time
dvdt_throat(k) = dvdt(throat);
drhodt_throat(k) = drhodt(throat);
%plotting the numerical results
%non-dimensional mass flow rate at the different time steps
figure(1)
if k == 1
plot(x,mass_flow_initial,'color','K','LineWidth', 1.3, 'LineStyle','--');
hold on
elseif k == 50
plot(x,mass_flow, 'm', 'LineWidth', 1.3);
hold on
elseif k == 100
plot(x,mass_flow, 'c', 'LineWidth', 1.3);
hold on
elseif k == 150
plot(x,mass_flow, 'g', 'LineWidth', 1.3);
hold on
elseif k == 200
plot(x,mass_flow, 'y', 'LineWidth', 1.3);
hold on
elseif k == 700
plot(x,mass_flow, 'k', 'LineWidth', 1.3);
hold on
end
title("Mass Flow Rate Distribution across the Nozzle [Non-Conservative Form]")
xlabel("Nozzle X-Distance")
ylabel("Non-Dimensional Mass Flow Rate")
legend('0^t^h Timestep', '50^t^h Timestep','100^t^h Timestep','150^t^h Timestep','200^t^h Timestep','700^t^h Timestep')
axis([0 3 0.1 2])
grid minor;
end
grid minor;
%plotting the density & velocity time derivative variation with time at the nozzle throat
figure(2)
plot(linspace(1,nt,nt),dvdt_throat, 'k', 'LineWidth', 1.3);
hold on
plot(linspace(1,nt,nt),drhodt_throat,'g', 'LineWidth', 1.3);
title("Time Derivatives of Velocity & Density [Non-Conservative Form]")
xlabel("Number of Time Steps")
ylabel("Residuals")
legend({'$dv/dt$' , '$drho/dt$'},'Interpreter', 'latex');
grid minor;
legend({'$dv/dt$' , '$partial{rho}/partial{t}$'},'Interpreter', 'latex');
%plotting steady-state Mach No. & Density along with the Nozzle
figure(3)
plot(x,rho,'k', 'LineWidth', 1.3);
xlabel('Nozzle X - Direction')
ylabel('Non-Dimensional Density')
yyaxis right
plot(x,Mach_no,'g', 'LineWidth', 1.3);
ylabel('Non-Dimensional Mach No.')
grid minor;
title("Steady State Density & Mach No. across Nozzle [Non-Conservative Form]")
legend('Density','Mach Number');
%printing total flow time
sprintf("Total Flow time for non conservative form = %.2f seconds", total_sim_time)
endCode for Conservative Form of ID Quasi Super-Sonic Nozzle Flow
function [v,pressure,T,Mach_no,mass_flow,rho,rho_throat,P_throat,v_throat,T_throat,Mach_no_throat, mass_flow_throat] = conservative_form(n,x,dx,gamma,CFL,nt);
%initial profile
%all flow variables are in their non-dimensional form
for i = 1:n
if 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);
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);
P(i) = rho(i)*T(i);
end
A(i) = 1 + 2.2*(x(i) - 1.5).^2; %converging diverging area
v(i) = 0.59/(rho(i)*A(i)); %velocity
mass_flow(i) = rho(i).*A(i).*v(i); %mass flow rate
end
mass_flow_initial = rho.*A.*v; %initial mass flow rate
throat = find(A == 1); %find throat location
%defining inital solution vector
U1 = rho.*A;
U2 = rho.*A.*v;
U3 = rho.*A.*((T./(gamma - 1)) + ((gamma/2).*v.^2));
total_sim_time = 0; %initial flow time
%Solving the equations (Conservative Form)
%Maccormack Method
for z = 1:nt; %time marching
% Storing the initial data
U1_old = U1;
U2_old = U2;
U3_old = U3;
%calculating the time step
dt = min((CFL.*dx)./(sqrt(T) + v));
%calculating the intial flux vector
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));
%PRIDECTOR STEP
for j = 2:n-1
%source term
J(j) = (1/gamma)*rho(j)*T(j)*((A(j+1) - A(j))/dx);
%continuity equation
dU1_dt_p(j) = -((F1(j+1) - F1(j))/dx);
%momemtum equation
dU2_dt_p(j) = -((F2(j+1) - F2(j))/dx) + J(j);
%energy equation
dU3_dt_p(j) = -((F3(j+1) - F3(j))/dx);
%updating the solution vector
U1(j) = U1(j) + dU1_dt_p(j)*dt;
U2(j) = U2(j) + dU2_dt_p(j)*dt;
U3(j) = U3(j) + dU3_dt_p(j)*dt;
end
%correcting flux vector
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));
%calculating premitive variables
rho = U1./A;
v = U2./U1;
T = (gamma - 1)*((U3./U1) - ((gamma/2)*(v.^2)));
P = rho.*T;
%CORRECTOR STEP
for k = 2:n-1
J(k) = (1/gamma)*rho(k)*T(k)*((A(k) - A(k-1))/dx); %source term
%continuity equation
dU1_dt_c(k) = -((F1(k) - F1(k-1))/dx);
%momemtum equation
dU2_dt_c(k) = -((F2(k) - F2(k-1))/dx) + J(k);
%energy equation
dU3_dt_c(k) = -((F3(k) - F3(k-1))/dx);
end
%calculating the averaged time derivatives
dU1_dt = 0.5*(dU1_dt_p + dU1_dt_c);
dU2_dt = 0.5*(dU2_dt_p + dU2_dt_c);
dU3_dt = 0.5*(dU3_dt_p + dU3_dt_c);
%calculating corrected solution vector
for m = 2:n-1
U1(m) = U1_old(m) + dU1_dt(m)*dt;
U2(m) = U2_old(m) + dU2_dt(m)*dt;
U3(m) = U3_old(m) + dU3_dt(m)*dt;
end
%apply BC
%inlet
U1(1) = rho(1)*A(1); %fixed value
U2(1) = 2*U2(2) - U2(3); %extrapolated
U3(1) = U1(1)*((T(1)/(gamma - 1)) + ((gamma/2)*(v(1)^2)));
%outlet
%extrapolated from the inside domain
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 corrected primitive non-dimensional variable for next time step
rho = U1./A;
v = U2./U1;
T = (gamma - 1)*((U3./U1) - ((gamma/2)*(v.^2)));
pressure = rho.*T;
mass_flow = rho.*A.*v;
Mach_no = (v./sqrt(T));
%calculating non-dimensioanl variables at the throat
rho_throat(z) = rho(throat);
v_throat(z) = v(throat);
T_throat(z) = T(throat);
P_throat(z) = P(throat);
mass_flow_throat(z) = mass_flow(throat);
Mach_no_throat(z) = Mach_no(throat);
total_sim_time = total_sim_time + dt; %updating the total flow time
%plotting the numerical results
%non-dimensional mass flow rate at the different time steps
figure(6)
if z == 1
plot(x,mass_flow_initial,'color','K','LineWidth', 1.3, 'LineStyle','--');
hold on
elseif z == 50
plot(x,mass_flow, 'm', 'LineWidth', 1.3);
hold on
elseif z == 100
plot(x,mass_flow, 'c', 'LineWidth', 1.3);
hold on
elseif z == 200
plot(x,mass_flow, 'g', 'LineWidth', 1.3);
hold on
elseif z == 400
plot(x,mass_flow, 'y', 'LineWidth', 1.3);
hold on
elseif z == 700
plot(x,mass_flow, 'k', 'LineWidth', 1.3);
hold on
end
title("Mass Flow Rate Distribution across the Nozzle [Conservative Form]")
xlabel("Nozzle X-Distance")
ylabel("Non-Dimensional Mass Flow Rate")
legend('0^t^h Timestep','50^t^h Timestep','100^t^h Timestep','200^t^h Timestep','400^t^h Timestep','700^t^h Timestep')
axis([0 3 0.4 1])
grid minor;
end
grid minor;
%printing total flow time
sprintf("Total Flow time for Conservative form = %.2f seconds", total_sim_time)
endResults & Discussion
The above graph shows the instantaneous non-dimensional mass flow rate distribution across the non-dimensional distance of the nozzle at the different time steps obtained by solving the non-conservative form of equations. There are 6 curves, each at different times during the time marching solution. The strange-looking dashed curve refers to the initial distribution(initial condition) of mass flow rate at time 0 second. After 50 time steps, the mass flow rate distribution inside the nozzle has changed considerably. After 100 time steps, the mass flow rate distribution inside the nozzle has changed radically, this is due to the transient variation of the flow field variable inside the nozzle. However, after 200 time-steps mass flow rate distribution begins to settle down, and after 700 time-steps mass flow rate distribution across the nozzle has become constant. This shows the mass flow has converged to steady-state value throughout the nozzle.
The above graph shows the residuals with respect to time obtained by solving the non-conservative form of equations. The time derivative of density and velocity are present in LHS of continuity and momentum equation, As time progresses and as steady-state is approached, the RHS of continuity and momentum equation should approach zero(derivatives of velocity and density should approach zero or velocity and density should not change with respect to time after steady-state) Since it is almost impossible to have right-hand-side precisely zero due to the nature of Numerical method, the y-axis is labeled as residuals.
At an early time derivatives of density and velocity are large, They oscillate in value. This oscillation is associated with various compression and expansion waves propagating through the nozzle during the transient process.
After a large time, around 500 time-steps, time derivatives settle down, become very small in the order of magnitude. This shows that flow variables are slowly reaching a steady-state value. Ideally, time derivatives of density and velocity should go to zero. However, this will never happen over a finite number of time steps.
The above graph shows the non-dimensional Mach no. & non-dimensional Density variation across the non-dimensional distance of the nozzle obtained by solving the non-conservative form of equations. We can see the flow is subsonic in the nozzle converging section(x<1.6), the flow becomes sonic at the nozzle throat(x = 1.6), Flow goes supersonic in diverging section of the nozzle(x > 1.6). Density decreases as flow expand through the nozzle.
The above graph shows the non-dimensional Density, Temperature, Pressure & Mach No. transient variation at the nozzle's throat location obtained by solving the non-conservative form of equations. The largest change takes place at early times, the transient behavior can be seen. The magnitude of the initial change in variable depends on applied initial conditions. But after a long course of time, around ~ 400 time steps, change in the variable is getting relaxed, and a steady-state is achieved.
The above graph shows the steady-state field of non-dimensional Density, Temperature, Pressure & Mach No. across the direction of the nozzles obtained by solving the non-conservative form of equations. The results show good agreement with the analytical results qualitatively and quantitatively(analytical results are not shown in the graph but can be easily obtained).
The above graph shows the instantaneous non-dimensional mass flow rate distribution across the non-dimensional distance of the nozzle at the different time steps obtained by solving the non0conservative form of equations. There are 6 curves, each at different times during the time marching solution. The dashed straight curve refers to the initial distribution(initial condition) of mass flow rate at time 0 second. The constant mass-flow rate = 0.59(close to analytical result) across the nozzle is assumed for the calculation of the non-conservative form of the equation. After the 700th time step, we can see the steady-state distribution of the mass-flow rate is achieved across the nozzle.
The below graph shows the steady-state mass flow rate profile(at a magnified scale) comparison between conservative and non-conservative forms numerical results after the 1000th time step when the convergence is achieved.
The dotted straight line shows the analytical results. The conservative form results give a distribution that is much closer to being constant and close to the analytical results. In contrast on a magnified scale, the Non-conservative form results show a sizable variation, with some spurious oscillation near the inflow and outflow. Of course, when plotted on a large scale these variation/oscillations are not apparent and the mass flow rate appears to be constant.
This shows the general advantage of the conservative form over the non-conservative form, the conservative form of the equation does a better job in preserving the mass flow rate throughout the nozzle, this is mainly because the mass-flow rate is the primary output from the conservative form of momentum equation which is not the case in the non-conservative form of the momentum equation.
Let's compare the Temperature, Density, Mach no. obtained analytically, from the conservative and non-conservative forms at the nozzle throat location.
Non-Conservative form numerical results are distinctly closer to the analytical results as compared with numerical results obtained from the conservative form of equations.
For the non-conservative form of numerical results, the mesh independence is achieved at 31 grid points itself which is not the case in numerical results obtained from the conservative form of equations.
After doubling the grid points numerical results obtained from the conservative form of equations are closer to the analytical results but still not as close as the non-conservative form of numerical results with half as many grid points.
The above graph shows the non-dimensional Density, Temperature, Pressure & Mach No. transient variation at the nozzle's throat location obtained by solving the conservative form of equations. If we compare this with numerical results obtained from non-conservative form, Here we can notice initial gradients/changes are not that too high and the flow variable started to relax after ~200 time steps whereas in the non-conservative form we can notice flow variable relaxes around ~400 time steps. Here faster convergence can be seen because of using better initial conditions.
The above graph shows the steady-state field of the Density, Pressure, Temperature & Mach no. across the nozzle X-Direction obtained by solving conservative form of governing equation. The results show good agreement with the analytical results qualitatively and quantitatively(analytical results are not shown in the graph but can be easily obtained).