Simulation of a 1D Super-sonic nozzle flow using Macormack Method (MATLAB)
Simulation of a 1D Super-sonic nozzle flow using Macormack Method
Introduction: In this project, we are going to simulate the isentropic flow through a quasi 1D subsonic-supersonic nozzle. We will derive both the conservation and non-conservation forms of the governing equations and solve them using the MacCormack's technique. Then we are going to determine the steady-state temperature distribution for the flow-field variables and investigate the difference between the two forms of governing equations by comparing their solutions.
The flow at the inlet to the nozzle comes from a reservoir where the pressure and temperature are denoted by `P_0`, and `T_0`, respectively. The cross-sectional area of the reservoir is large, and hence the velocity is very small (V ~ 0). Thus, `P_0` and `T_0` are the stagnation values, or total pressure and total temperature, respectively. The flow expands isentropically to supersonic speeds at the nozzle exit, where the exit pressure, temperature, velocity, and Mach number are denoted by `P_e, T_e,V_e` and `M_e` respectively. We assume that flow properties vary only along the x-direction. Such flow is defined as quasi-one-dimensional flow.
Derivation of Governing Equations:
Continuity Equation- The integral form of the governing equation can be given as, `del/(delt)int int int_V rhodV + int int_S rho V.dS = 0` where, `int int_Srho V.dS` is equal to Outflow - Inflow. We will apply this equation to the shaded control volume.
The volume integral becomes, `del/(delt)int int int_V rho. dV = del/(delt)int int int rho A.dx = del/(delt)rho A.dx`
The surface integral becomes, `int int_S rho V.dS = - rhoAV + (rho+drho)(V+dV)(A+dA)`
`= -rhoAV + rhoVA + rhoVdA +rhoAdV + rhodVdA + AVdrho + VdAdrho +VdAdrho + AdVdrho+drhodVdA`
neglecting the product of diffrentials since they are very small, our equation becomes, `intint_SrhoV.dS = rho VdA + rhoAdV +AVdrho = d(rhoAV)`
Substituting these volume and surface integral in continuity equation we get, `del/(delt)(rhoAdx) + d(rhoAV) = 0`
Dividing the above equation by dx, `(del(rhoA))/(delt) + (del(rhoAV))/(delx) = 0` this is quazi-one dimensional continuity equation.
Non-Conservative form continuity equation:
Expanding the differential, `(del(rhoA))/(delt) +V(del(rhoA))/(delx) + rhoA(delV)/(delx) = 0`
Momentum Equation -
Integral form of momentum equation is given by, `del/(delt)intintint_Vrhou.dV + intint_S(rho u V).dS = -intint_S(pdS)_x`
Applying this to our control volume as we have done for the continuity equation,
The volume integral becomes, `del/(delt)intintint_Vrhou.dV = del/(delt)(rhoVAdx)`
The surface integral becomes, `intint_S(rho u V).dS = - rhoAV^2 + (rho+drho)(V+dV)^2(A+dA)`
The pressure force term, `int int (pdS)_x = -pA + (p+dp)(A+dA)- 2p((dA)/2)`
Substituting these terms in momentum equation, `del/(delt)(rhoVAdx) - rhoV^2A + (rho+drho)(V+dV)^2(A+dA) = pA - (p+dp)(A+dA) + pdA`
neglecting the product of diffrentials since they are very small, our equation becomes,
`del/(delt)(rhoVAdx) + d(rhoV^2A) = -Adp`
Dividing the above equation by dx, `(del(rhoVA))/(delt) + (del(rhoV^2A))/(delx) = -A(delp)/(delx)` this is the conservative form of momentum equation.
Non-Conservative form:
We will multiply the continuity equation by V, `V*(del(rhoA))/(delt) + V*(del(rhoAV))/(delx) = 0` and we will subtract this equation from conservative form of momentum equation.
`(del(rhoVA))/(delt) + (del(rhoV^2A))/(delx) - V*(del(rhoA))/(delt) - V*(del(rhoAV))/(delx) = -A(delp)/(delx)`
`V(del(rhoA))/(delt) + Arho(del(V))/(delt) + V^2(del(rhoA))/(delx) + 2rhoA(del(V))/(delx) - V*(del(rhoA))/(delt) - VrhoA(del(V))/(delx) - V^2(del(rhoA))/(delx)= -A(delp)/(delx)`
Cancelling the terms we get, `rho(delV)/(delt) + rhoV(delV)/(delx) = -(delp)/(delx)`. This is the Non-Conservative form of momentum equation.
Energy Equation -
Integral form of energy equation, `del/(delt)intintint_V rho(e+ (V^2)/2)dV + intint_S rho (e + V^2/2)V.dS = - intint_SrhoV.dS`
Allpying this to the control volume as we have done for the above equations and neglecting the product of differentials, we get the conservative form of energy equation.
`(del[rho(e+ V^2/2)A])/(delt) + (del[rho(e+V^2/2)A])/(delx) = - (del(pAV))/(delx)`
Non-coservative form of energy equation:
The non-coservative equation is obtained by multiplying continuity equation by 'e' and subtracting it from the conservative form of the energy equation as we have done for momentum equation.
`rho(dele)/(delt) + +rhoV(dele)/(delx) = -rho(delV)/(delx) - rhoV(del(lnA))/(delx)`
Non-Dimensional Form of Governing Equations:
We will replace `p =rhoRT and (delp)/(delx) = R(rho(delT)/(delx) + T(delrho)/(delx))`. `e = c_v*T`
The non-dimensional temperature, density, velocity, length, time , area is `T^' = T/T_0 , rho^' = rho/rho_0, V^' = V/a_0, x^' = x/L. t^' = t/(L/a_0), A^' = A/(A^+)`
Where `T_0 and rho_0` are properties at reservoir. `A^+` is the sonic throat, `a_0` is speed of sound and `L` is length of nozzel.
Conservative form:
Continuity Equation- `(del(rho^' A^'))/(delt^') + (del(rho^'A^'V^'))/(delx^' ) = 0`
Momentum Equation- `(del(rho^'A^'V^'))/(delt^') + (del(rho^'A^'V^('2) + (1/gamma)p^'A^'))/(dx) = 1/gammap^'(delA^')/(delx)`
Energy Equation- `(del[rho^'(e^'/(gamma-1) + gamma/2 V^('2))A^'])/(dt^') + (del[rho^'(e^'/(gamma-1) + gamma/2 V^('2))A^'V^' + rho^'A^'V^' ])/(dx^') = 0`
Non Conservative form:
Continuity Equation- `(delrho^')/(delt) = -rho^'(delV^')/(delx^')-rho^' V^'(delln(A^'))/(delx^') - V^'(delrho^')/(delx^')`
Momentum equation- `(delV^')/(delt^') = -V^'(delV^')/(delx^') - 1/gamma((delT^')/(delx^') + T^'/rho^'(delrho^')/(delx^'))`
Energy Equation- `(delT^')/(delt^') = -V^'(delT^')/(delx^')-(gamma-1)T^'[(delV^')/(delx^')+V^'(del(lnA^'))/(delx^')]`
Nozzle Shape: `A = 1+2.2(x-1.5)^2` `0lexle3`
Initial Condition at t = 0:
Non-Conservative - `rho = 1-0.3146x` `T = 1-0.2314x` `V = (0.1 + 1.09x)T^(1/2)`
Conservative -
for `0 le x le 0.5`, `rho` = 1; T = 1
for `0.6 le x le 1.5` `rho` = 1 - 0.366*(x - 0.5); T = 1 - 0.167*(x - 0.5)
for `1.6lexle3` `rho` = 0.634 - 0.3879*(x - 1.5); T = 0.833 - 0.3507*(x - 1.5)
`V = 0.59/(rho*A)`
MacCormack Method:
MacCormack method is a discretization scheme for the numerical solution of hyperbolic partial differential equations. This is a second-order finite difference method that works for non-linear hyperbolic problems. The CFL criteria (explicit scheme) is `sigma le 1` for linear equation. The MacCormack method is a two-step method, the first step is the Predictor step and second is the Corrector step. It is a variant of the Lax-Wendroff approach.
From the Von-Neumann analysis for non-linear equation (Linearize) we get `stability (Deltat)/(Deltax) < 1/(|u|`
Let us consider a convection equation, `(delrho)/(delt) =-(rho(delu)/(delx) +u(delrho)/(delx))`, our aim is to calculate the value of `rho` at `t+Deltat`.
This is given by, `rho_i^(t+Deltat) = rho_i^t + ((delrho)/(delt))_(avg)Deltat`, where `((delrho)/(delt))_(avg)` is the mean representative value between `t` and `t+Deltat`. `((delrho)/(delt))_(avg)` is calculated so as to preserve the second order accuracy without the need to calculate the second order time derivative.
Step_1 (Predictor, Forward difference in space): `((delrho)/(delt))_i^t = -(rho_i(u_(i+1) - u_i)/(Deltax) + u_i (rho_(i+1) - rho_i)/(Deltax))^t`, from this step the predicted value of `barrho_i^(t+Deltat)` is calculated using Taylor series (first two terms). This value is first order accurate.
`barrho_i^(t+Deltat) = rho_i^t + ((delrho)/(rhot))_i^tDeltat`
Step_2 (Corrector, Backward difference in space with `(Deltat)/2`): `((delbarrho)/(delt))_i^(t+Deltat) = -(barrho_i(baru_i - baru_(i-1))/(Deltax) + baru_i (barrho_i - barrho_(i-1))/(Deltax))^(t+Deltat)`.
Average time derivative `((delrho)/(delt))_(avg) = 1/2 [ ((delrho)/(delt))_i^t + ((delbarrho)/(delt))_i^(t+Deltat)]`
`rho_i^(t+Deltat) = rho_i^t + ((delrho)/(delt))_(avg)Deltat`
We are going to apply this method to our governing equations.
Program: Non-Conservative
Function: non_conservative
function [v, T,rho] = non_conservative()
global N_Mass_Flow % non-conservative; mass flow rate at throat
z = 1;
n_iter = [31 61];
for u = 1:length(n_iter)
n = n_iter(u);
x_throat = (n+1)/2;
% input
x = linspace(0,3,n);
dx = x(2) - x(1);
gamma = 1.4;
% calculate initial conditions at t = 0
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 = A/A^*
a = (T).^0.5; % speed of sound (non-dimensional)
% time step
nt = 1400;
C = 0.5;
dt = min((C*dx)./(a+ v));
rho_old = rho;
v_old = v;
T_old = T;
time_steps = linspace(1,nt,nt);
% outer time loop
for k = 1:nt
%---------[ Predictor step ]---------%
for i = 2:n-1
dvdx(i) = (v(i+1) - v(i))/dx;
drhodx(i) = (rho(i+1) - rho(i))/dx;
dlnAdx(i) = (log(A(i+1)) - log(A(i)))/dx;
dTdx(i) = (T(i+1) - T(i))/dx;
% Continuity equation
drhodt_p(i) = -rho(i)*dvdx(i) - rho(i)*v(i)*dlnAdx(i) - v(i)*drhodx(i);
% Momentum equation
dvdt_p(i) = -v(i)*dvdx(i) - (1/gamma)*(dTdx(i) + (T(i)/rho(i))*drhodx(i));
% Energy equation
dTdt_p(i) = -v(i)*dTdx(i) - (gamma - 1)*T(i)*(dvdx(i) + v(i)*dlnAdx(i));
% Solution update
[rho(i)] = rho(i) + drhodt_p(i)*dt;
v(i) = v(i) + dvdt_p(i)*dt;
T(i) = T(i) + dTdt_p(i)*dt;
end
%---------[ Corrector Step ]--------%
for j = 2:n-1
dvdx(j) = (v(j) - v(j-1))/dx;
drhodx(j) = (rho(j) - rho(j-1))/dx;
dlnAdx(j) = (log(A(j)) - log(A(j-1)))/dx;
dTdx(j) = (T(j) - T(j-1))/dx;
% Continuity equation
drhodt_c(j) = -rho(j)*dvdx(j) - rho(j)*v(j)*dlnAdx(j) - v(j)*drhodx(j);
% Momentum equation
dvdt_c(j) = -v(j)*dvdx(j) - (1/gamma)*(dTdx(j) + (T(j)/rho(j))*drhodx(j));
% Energy equation
dTdt_c(j) = -v(j)*dTdx(j) - (gamma - 1)*T(j)*(dvdx(j) + v(j)*dlnAdx(j));
end
% Average Solution
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 q = 2:n-1
rho(q) = rho_old(q) + drhodt(q)*dt;
v(q) = v_old(q) + dvdt(q)*dt;
T(q) = T_old(q) + dTdt(q)*dt;
end
% Apply boundary condition
v(1) = 2*v(2) - v(3); % Inlet
v(n) = 2*v(n-1) - v(n-2); % Outlet
rho(n) = 2*rho(n-1) - rho(n-2); % Outlet
T(n) = 2*T(n-1) - T(n-2); % Outlet
% Value(Non-D) at throat at diff time step, A = 1
throat_T(k) = T(x_throat)./ T(1);
throat_v(k) = v(x_throat)./((T(1)).^0.5);
throat_rho(k) = rho(x_throat)./rho(1);
% Update
rho_old = rho;
v_old = v;
T_old = T;
dt = min((C*dx)./(a+ v));
if k == 100 || k==200 || k==400 || k==700|| k==1000
Mass_Flow = (rho.*A.*v);
N_Mass_Flow{z} = Mass_Flow./(rho(1)*A(16)*((T(1)).^0.5));
z = z + 1;
end
end
if n == 31
figure(1)
plot(x, N_Mass_Flow{1},'--','LineWidth',1.3)
hold on
plot(x, N_Mass_Flow{2},'--','LineWidth',1.3)
plot(x, N_Mass_Flow{3},'--','LineWidth',1.3)
plot(x, N_Mass_Flow{4},'--','LineWidth',1.3)
plot(x, N_Mass_Flow{5},'--','LineWidth',1.3)
grid on
title('NC* Mass Flow at Different Time Step', 'FontSize', 12)
legend('100dt','200dt','400dt','700dt','1000dt');
xlabel('Nondimensional Distance Through Nozzel', 'FontSize', 12)
ylabel('Mass Flow Rate','FontSize', 12)
figure(2)
plot(time_steps, throat_T, 'r', 'LineWidth',1.3)
hold on
grid on
plot(time_steps, throat_v, 'g--', 'LineWidth',1.3)
plot(time_steps, throat_rho, 'b.', 'LineWidth',1.3)
title('NC* Variables at Nozzel Throat, A = 1','FontSize', 12)
legend('Temperature', 'Velocity', 'Density')
xlabel('Time Steps', 'FontSize', 12)
ylabel('Variables (Non-Dimensional)','FontSize', 12)
figure(3)
plot(x,v,'g--', 'LineWidth',1.3)
hold on
grid on
plot(x, T, 'r', 'LineWidth',1.3)
plot(x, rho, 'b.', 'LineWidth',1.3)
title('NC* Steady State Distribution Inside Nozzel', 'FontSize', 12)
legend('Velocity', 'Temperature', 'Density')
xlabel('Nondimensional Distance Through Nozzel', 'FontSize', 12)
ylabel('Variables (Non-Dimensional)','FontSize', 12)
end
Density = "(Non-Conservative) Density for " +n +" points = " +throat_rho(1000)
Temperature = "(Non-Conservative) Temperature for " +n +" points = " +throat_T(1000)
Velocity = "(Non-Conservative) Velocity for " +n +" points = " +throat_v(1000)
end
end
Function: conservative
function [v, T,rho] = conservative()
global C_Mass_Flow % conservative; mass flow rate at throat
z = 1;
n_iter = [31 61];
for u = 1:length(n_iter)
n = n_iter(u);
x_throat = (n+1)/2;
% input
x = linspace(0,3,n);
dx = x(2) - x(1);
gamma = 1.4;
%----calculate initial conditions at t = 0----%
x1 = (0.5/dx) + 1;
for xx = 1:x1
rho(xx) = 1; % Density
T(xx) = 1; % Temperature
end
x2 = (1.5/dx) + 1;
for xx = x1:x2
rho(xx) = 1 - 0.366*(x(xx) - 0.5);
T(xx) = 1 - 0.167*(x(xx) - 0.5);
end
x3 = (3/dx) + 1;
for xx = x2:x3
rho(xx) = 0.634 - 0.3879*(x(xx) - 1.5);
T(xx) = 0.833 - 0.3507*(x(xx) - 1.5);
end
e = T; % Energy
A = 1 + 2.2*(x-1.5).^2; % area = A/A^*
v = 0.59./(rho.*A); % Velocity
a = (T).^0.5; % speed of sound (non-dimensional)
p = rho.*T; % Pressure
p_old = p;
%--------------------------------------------%
% time step
nt = 1400;
C = 0.5;
dt = min((C*dx)./(a+ v));
%dt = 0.0001;
time_steps = linspace(1,nt,nt);
%----Solution vector U, Flux vector F, Source J----%
U1 = rho.*A;
U2 = 0.59.*ones(1,n);
U3 = rho.*((T/(gamma - 1)) + (gamma/2)*v.^2).*A; % e = T;
F1 = rho.*v.*A;
F2 = rho.*A.*v.^2 + ((p.*A)/gamma);
F3 = rho.*(T./(gamma-1) + (v.^2*gamma)/(2)).*v.*A + p.*A.*v; % e = T;
U1_old = U1;
U2_old = U2;
U3_old = U3;
F1_old = F1;
F2_old = F2;
F3_old = F3;
rho_old = rho;
v_old = v;
T_old = T;
%--------------------------------------------------%
% Outer Time Loop
for k = 1:nt
%---------[ 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;
dAdx(i) = (A(i+1) - A(i))/dx;
J2(i) = (p(i)*dAdx(i))/gamma;
% Continuity equation
dU1dt_p(i) = -dF1dx(i);
% Momentum equation
dU2dt_p(i) = -dF2dx(i) + J2(i);
% Energy equation
dU3dt_p(i) = -dF3dx(i);
% Solution update
U1(i) = U1(i) + dU1dt_p(i)*dt;
rho(i) = U1(i)/A(i);
U2(i) = U2(i) + dU2dt_p(i)*dt;
v(i) = U2(i)/U1(i);
U3(i) = U3(i) + dU3dt_p(i)*dt;
T(i) = (gamma - 1)*((U3(i)/U1(i)) - (gamma/2)*v(i)^2);
p(i) = rho(i)*T(i);
F1(i) = rho(i)*v(i)*A(i);
F2(i) = rho(i)*A(i)*v(i)^2 + ((p(i)*A(i))/gamma);
F3(i) = rho(i)*(T(i)/(gamma-1) + (v(i)^2*gamma)/(2))*v(i)*A(i) + p(i)*A(i)*v(i);
end
%---------[ Corrector Step ]--------%
for j = 2:n-1
dF1dx(j) = (F1(j) - F1(j-1))/dx;
dF2dx(j) = (F2(j) - F2(j-1))/dx;
dF3dx(j) = (F3(j) - F3(j-1))/dx;
dAdx(j) = (A(j) - A(j-1))/dx;
J2(j) = (p(j)*dAdx(j))/gamma;
% Continuity equation
dU1dt_c(j) = -dF1dx(j);
% Momentum equation
dU2dt_c(j) = -dF2dx(j) + J2(j);
% Energy equation
dU3dt_c(j) = -dF3dx(j);
end
% Average Solution
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 q = 2:n-1
U1(q) = U1_old(q) + dU1dt(q)*dt;
rho(q) = U1(q)/A(q);
U2(q) = U2_old(q) + dU2dt(q)*dt;
v(q) = U2(q)/U1(q);
U3(q) = U3_old(q) + dU3dt(q)*dt;
T(q) = (gamma - 1)*((U3(q)/U1(q)) - (gamma/2)*v(q)^2);
p(q) = rho(q)*T(q);
F1(q) = rho(q)*v(q)*A(q);
F2(q) = rho(q)*A(q)*v(q)^2 + ((p(q)*A(q))/gamma);
F3(q) = rho(q)*(T(q)/(gamma-1) + (v(q)^2*gamma)/(2))*v(q)*A(q) + p(q)*A(q)*v(q);
a(q) = T(q)^0.5;
end
% Apply boundary condition
U1(1) = rho(1)*A(1);
U2(1) = 2*U2(2) - U2(3);
v(1) = U2(1)/U1(1);
U3(1) = rho(1)*((T(1)/(gamma - 1)) + (gamma/2)*v(1)^2)*A(1);
F1(1) = rho(1)*v(1)*A(1);
F2(1) = rho(1)*A(1)*v(1)^2 + ((p(1)*A(1))/gamma);
F3(1) = rho(1)*(T(1)/(gamma-1) + (v(1)^2*gamma)/(2))*v(1)*A(1) + p(1)*A(1)*v(1);
U1(n) = 2*U1(n-1) - U1(n-2);
rho(n) = U1(n)/A(n);
U2(n) = 2*U2(n-1) - U2(n-2);
v(n) = U2(n)/U1(n);
U3(n) = 2*U3(n-1) - U3(n-2);
T(n) = (gamma - 1)*((U3(n)/U1(n)) - (gamma/2)*v(n)^2);
p(n) = rho(n)*T(n);
F1(n) = rho(n)*v(n)*A(n);
F2(n) = rho(n)*A(n)*v(n)^2 + ((p(n)*A(n))/gamma);
F3(n) = rho(n)*(T(n)/(gamma-1) + (v(n)^2*gamma)/(2))*v(n)*A(n) + p(n)*A(n)*v(n);
a(n) = T(n)^0.5;
% Value(Non-D) at throat at diff time step, A = 1
throat_T(k) = T(x_throat)./ T(1);
throat_v(k) = v(x_throat)./((T(1)).^0.5);
throat_rho(k) = rho(x_throat)./rho(1);
% Update
U1_old = U1;
U2_old = U2;
U3_old = U3;
dt = min((C*dx)./(a+ v));
if k == 100 || k==200 || k==400 || k==700|| k==1400
Mass_Flow = (rho.*A.*v);
C_Mass_Flow{z} = Mass_Flow./(rho(1)*A(16)*((T(1)).^0.5));
z = z + 1;
end
end
if n == 31
figure(4)
plot(x, C_Mass_Flow{1},'--','LineWidth',1.3)
hold on
plot(x, C_Mass_Flow{2},'--','LineWidth',1.3)
plot(x, C_Mass_Flow{3},'--','LineWidth',1.3)
plot(x, C_Mass_Flow{4},'--','LineWidth',1.3)
plot(x, C_Mass_Flow{5},'--','LineWidth',1.3)
grid on
title('C* Mass Flow at Different Time Step', 'FontSize', 12)
legend('100dt','200dt','400dt','700dt','1000dt');
xlabel('Nondimensional Distance Through Nozzel', 'FontSize', 12)
ylabel('Mass Flow Rate','FontSize', 12)
figure(5)
plot(time_steps, throat_T, 'r', 'LineWidth',1.3)
hold on
grid on
plot(time_steps, throat_v, 'g--', 'LineWidth',1.3)
plot(time_steps, throat_rho, 'b.', 'LineWidth',1.3)
title('C* Variables at Nozzel Throat, A = 1','FontSize', 12)
legend('Temperature', 'Velocity', 'Density')
xlabel('Time Steps', 'FontSize', 12)
ylabel('Variables (Non-Dimensional)','FontSize', 12)
figure(6)
plot(x,v,'g--', 'LineWidth',1.3)
hold on
grid on
plot(x, T, 'r', 'LineWidth',1.3)
plot(x, rho, 'b.', 'LineWidth',1.3)
title('C* Steady State Distribution Inside Nozzel', 'FontSize', 12)
legend('Velocity', 'Temperature', 'Density')
xlabel('Nondimensional Distance Through Nozzel', 'FontSize', 12)
ylabel('Variables (Non-Dimensional)','FontSize', 12)
end
Density = "(Conservative) Density for " +n +" points = " +throat_rho(1000)
Temperature = "(Conservative) Temperature for " +n +" points = " +throat_T(1000)
Velocity = "(Conservative) Velocity for " +n +" points = " +throat_v(1000)
end
end
Main Program:
clear all ; close all; clc
% [ 1D SUPER-SONIC NOZZLE FLOW SIMULATION USING MACORMACK METHOD ]
global C_Mass_Flow N_Mass_Flow
n = 31;
x = linspace(0,3,n);
Analytical = 0.579*ones(n);
Non_Conservative = non_conservative();
Conservative = conservative();
figure(7)
plot(x,N_Mass_Flow{5}, 'LineWidth',1.3)
hold on
plot(x,C_Mass_Flow{5},'LineWidth',1.3)
plot(x,Analytical,'k--', 'LineWidth',1.3)
grid on
title('Mass Flow Rate Distribution Inside Nozzle', 'FontSize', 12)
legend('Non-Conservative','Conservative','Analytical');
xlabel('Nondimensional Distance Through Nozzel', 'FontSize', 12)
ylabel('Mass Flow Rate','FontSize', 12)
Output:
-------------------------------------------------------------------------------------------Non-Conservative-------------------------------------------------------------------------------------------
Conditions at Nozzel Throat
| `rho^ast/rho_0` | `T^ast/T_0` | `M` | |
| Case 1: 31 Points | 0.63869 | 0.83646 | 0.91401 |
| Case 2: 61 Points | 0.6393 | 0.83625 | 0.91365 |
| Exact Analytical Solution | 0.634 | 0.833 | 1 |
----------------------------------------------------------------------------------------------Conservative----------------------------------------------------------------------------------------------
Conditions at Nozzel Throat
| `rho^ast/rho_0` | `T^ast/T_0` | `M` | |
| Case 1: 31 Points | 0.65059 | 0.84011 | 0.89948 |
| Case 2: 61 Points | 0.63857 | 0.83529 | 0.91329 |
| Exact Analytical Solution | 0.634 | 0.833 | 1 |
----------------------------------------------------------------------------------------------Comparision----------------------------------------------------------------------------------------------
Conclusion:
In this study, we have performed a numerical simulation of the 1-D supersonic nozzle using MacCormack method. The results obtained from a conservative form of governing equations are more accurate than a non-conservative form of governing equations as compared with the exact analytical solutions. The discretization of the non-conservative form of equations leads to 'internal sources'. The conservative and non-conservative form of equations are mathematically equivalent for arbitrary nonlinear fluxes but their numerical implementation is not equivalent!