Week 7 - Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
Aim: Simulating quasi 1-D supersonic flow through a converging-diverging nozzle using the MacCormack method.
Objectives:
In this challenge, you are going to simulate the isentropic flow through a quasi-1D subsonic-supersonic nozzle. You will derive both the conservation and non-conservation forms of the governing equations and solve them using MacCormack's technique. You need 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.
We need to show the following plots inside your report
1. Steady-state distribution of primitive variables inside the nozzle
2. Time-wise variation of the primitive variables
3. Variation of Mass flow rate distribution inside the nozzle at different time steps during the time-marching process
4. Comparison of Normalized mass flow rate distributions of both forms
Expected results:
- To perform a grid independence test on the solution and explain all findings.
- The code must be modular. You need to write separate functions to solve both forms.
- Discuss all the plots you are providing with observed inferences.
Theory:
A nozzle is a type of device designed to control the direction or characteristics of fluid flow at the exit of a pipe (or such a system). A nozzle is often a pipe or tube of varying cross-sectional area and thus can be used to direct or modify the flow of a fluid.
A convergent-divergent nozzle (also known as CD nozzle or de Laval nozzle) is a device used to accelerate a hot, pressurized gas passing through it to a higher supersonic speed in the axial direction by converting the heat energy of the flow to kinetic energy. Hence, such nozzles are widely used in steam turbines, supersonic jet engines, and rocket engine nozzles. Its operation relies on the different properties of gases flowing at subsonic, sonic, and supersonic speeds.
In the field of CFD, the MacCormack method is a prominent discretization scheme for obtaining numerical solutions to hyperbolic partial differential equations. (based on the second-order finite difference method). It is the variation of the two-step Lax-Wendroff scheme but much simple in application.
In this project, the problem is taken from JD Anderson's book on CFD.
The symmetrical area function for the CD nozzle is given as,
(AA∗)2=1+2.2(x−1.5)2(AA∗)2=1+2.2(x-1.5)2 @ 0≤x≤30≤x≤3
where A∗A∗ is the area at the throat of the nozzle and xx is the distance along the length of the nozzle.
The non-dimensionalized form of the variables required for this problem is,
x'=xL; V'=Va0; t'=tLa0; ρ'=ρρ0; T'=TT0; p'=pp0; A'=AA∗x′=xL; V′=Va0; t′=tLa0; ρ′=ρρ0; T′=TT0; p′=pp0; A′=AA∗
Governing Equations:
1. Non-Conservative form
The non-conservative form of the governing equation which represents continuity, momentum, and energy in terms of non-dimensionless variables,
Continuity: ∂ρ'∂t'=−ρ'∂V'∂x'−ρ'V'∂(lnA')∂x'−V'∂ρ'∂x'∂ρ′∂t′=-ρ′∂V′∂x′-ρ′V′∂(lnA′)∂x′-V′∂ρ′∂x′
Momentum: ∂V'∂t'=−V'∂V'∂x'−1γ[∂T'∂x'+T'ρ'∂ρ'∂x']∂V′∂t′=-V′∂V′∂x′-1γ[∂T′∂x′+T′ρ′∂ρ′∂x′]
Energy: ∂T'∂t'=−V'∂T'∂x'−(γ−1)T'[∂V'∂x'+V'∂(lnA')∂x']∂T′∂t′=-V′∂T′∂x′-(γ-1)T′[∂V′∂x′+V′∂(lnA′)∂x′]
The initial conditions as provided in the problem statement were,
ρ'(t=0)=1−0.3146x'ρ′(t=0)=1-0.3146x′
T'(t=0)=1−0.2314x'T′(t=0)=1-0.2314x′
V'(t=0)=(0.1+1.09x')√T'(t=0)V′(t=0)=(0.1+1.09x′)T′(t=0)
2. Conservative Form
The conservative form of governing equation which represents continuity, momentum, and energy in non-dimensionalized form is given as,
Continuity: ∂(ρ'A')∂t'+∂(ρ'A'V')∂x'=0∂(ρ′A′)∂t′+∂(ρ′A′V′)∂x′=0
Momentum: ∂(ρ'A'V')∂t'+∂(ρ'A'V'2+(1γ)p'A')∂x'=1γp'∂A'∂x'∂(ρ′A′V′)∂t′+∂(ρ′A′V′2+(1γ)p′A′)∂x′=1γp′∂A′∂x′
Energy: ∂[ρ'(T'γ−1+γ2V'2)A']∂t'+∂[ρ'(T'γ−1+γ2V'2)V'A'+p'A'V']∂x'=0∂[ρ′(T′γ-1+γ2V′2)A′]∂t′+∂[ρ′(T′γ-1+γ2V′2)V′A′+p′A′V′]∂x′=0
To simplify these equations, we substitute the following terms, which include the solution vectors UiUi, the flux vectors FiFi, and the source term J2J2.
U1=ρ'A'U1=ρ′A′
U2=ρ'A'V'U2=ρ′A′V′
U3=ρ'[T'γ−1+γ2V'2]A'U3=ρ′[T′γ-1+γ2V′2]A′
F1=ρ'A'V'=U2F1=ρ′A′V′=U2
F2=ρ'A'V'2+1γp'A'=U22U1+γ−1γ[U3−γ2U22U1]F2=ρ′A′V′2+1γp′A′=U22U1+γ-1γ[U3-γ2U22U1]
F3=ρ'(T'γ−1+γ2V'2)V'A'+p'A'V'=γU2U3U1−γ(γ−1)2U32U21F3=ρ′(T′γ-1+γ2V′2)V′A′+p′A′V′=γU2U3U1-γ(γ-1)2U23U12
J2=1γp'∂A'∂x'J2=1γp′∂A′∂x′
After substituting these terms in governing equations, we have,
∂U1∂t'=−∂F1∂x'∂U1∂t′=-∂F1∂x′
∂U2∂t'=−∂F2∂x'+J2∂U2∂t′=-∂F2∂x′+J2
∂U3∂t'=−∂F3∂x'∂U3∂t′=-∂F3∂x′
The flow variables are calculated on the numerical solution of the above terms,
ρ'=U1A'ρ′=U1A′
V'=U2U1V′=U2U1
T'=(γ−1)[U3U1−γ2V'2]T′=(γ-1)[U3U1-γ2V′2]
p'=ρ'T'p′=ρ′T′
The initial conditions as provided in the problem statement are,
ρ'(t=0)=⎧⎪⎨⎪⎩1.0for0≤x'<0.5(1.0−0.366(x'−0.5))for0.5≤x'<1.5(0.634−0.3879(x'−1.5))for1.5≤x'≤3ρ′(t=0)={1.0for0≤x′<0.5(1.0-0.366(x′-0.5))for0.5≤x′<1.5(0.634-0.3879(x′-1.5))for1.5≤x′≤3
T'(t=0)=⎧⎪⎨⎪⎩1.0for0≤x'<0.5(1.0−0.167(x'−0.5))for0.5≤x'<1.5(0.833−0.3507(x'−1.5))for1.5≤x'≤3T′(t=0)={1.0for0≤x′<0.5(1.0-0.167(x′-0.5))for0.5≤x′<1.5(0.833-0.3507(x′-1.5))for1.5≤x′≤3
V'(t=0)=U2ρ'(t=0)A'=0.59ρ'(t=0)A'V′(t=0)=U2ρ′(t=0)A′=0.59ρ′(t=0)A′ (considering a constant mass flow rate of U2=ρ'A'V'=0.59U2=ρ′A′V′=0.59 just for the initial condition)
MacCormack method of solving:
The MacCormack numerical scheme is a common method to solve flow equations. It consists of four steps:
- Predictor step
- Intermediate solution update
- Corrector step
- Final solution update
In the predictor step, the spatial derivative is calculated using the forward difference scheme. This time derivative from the governing equations obtained is recorded as ∂u(∂t)p∂u(∂t)p. Consider, the continuity equation for the non-conservative form. From the predictor step, we obtain,
(∂ρi∂t)P=−ρ'i[V'i+1−V'iΔx']−ρ'iV'i[(ln(A'i+1−ln(A'i))Δx']−V'i[ρ'i+1−ρ'iΔx'](∂ρi∂t)P=-ρ′i[V′i+1-V′iΔx′]-ρ′iV′i[(ln(A′i+1-ln(A′i))Δx′]-V′i[ρ′i+1-ρ′iΔx′]
During the intermediate solution update, the variable is updated based on the time derivative obtained from the predictor step as:
ρt+Δti,P=ρ'i+(∂ρ'i∂t')PΔtρi,Pt+Δt=ρ′i+(∂ρ′i∂t′)PΔt
Using these updated values, the governing equations are solved again using backward differencing in space in the corrector step. The time-derivative obtained during the corrector step is stored as (∂u∂t)C(∂u∂t)C.
(∂ρi∂t)C=−ρ'i,P[V'i,P−V'i−1,PΔx']−ρ'i,PV'i,P[(ln(A'i,P−ln(A'i−1,P))Δx']−V'i,P[ρ'i,P−ρ'i−1,PΔx'](∂ρi∂t)C=-ρ′i,P[V′i,P-V′i-1,PΔx′]-ρ′i,PV′i,P[(ln(A′i,P-ln(A′i-1,P))Δx′]-V′i,P[ρ′i,P-ρ′i-1,PΔx′]
For the final solution update, we use the average of the time-derivatives obtained using the predictor and the corrector steps. Since we are using the forward difference in the predictor step and backward difference in the corrector step, the average of the time-derivative for the final update, the overall scheme becomes a second-order numerical scheme.
(∂ρ∂t)av=12[(∂ρ∂t)P+(∂ρ∂t)C](∂ρ∂t)av=12[(∂ρ∂t)P+(∂ρ∂t)C]
ρt+Δti=ρti+(∂ρ∂t)avΔtρit+Δt=ρit+(∂ρ∂t)avΔt
Matlab Code:
The code is written in MATLAB.
1. Main Script
close all
clc
% Initializing the domain
L = 3;
nx = 31;
x = linspace(0,L,nx)';
dx = x(2)-x(1);
C = 0.5;
nt = 1400;
gamma = 1.4;
% Nozzle area profile
A = 1 + 2.2*(x-1.5).^2;
% Computing the flow variable through conservative and non-conservative form
[rho1, V1, T1, p1, M1, m1, rho_t1, V_t1, p_t1, T_t1, M_t1, residual1] = nonconservative(x, dx, A, C, nt, gamma);
[rho2, V2, T2, p2, M2, m2, rho_t2, V_t2, p_t2, T_t2, M_t2, residual2] = conservative(x, dx, A, C, nt, gamma);
% Plotting the results
figure(1)
subplot(4,1,1)
plot(x,M1)
grid on; ylim([0 4.5]); yticks(0:1.5:4.5);
ylabel('Mach Number')
title('Flow parameter (non-conservative form of governing equation)')
subplot(4,1,2)
plot(x,p1)
grid on; yticks(0:0.2:1);
ylabel('Pressure Ratio')
subplot(4,1,3)
plot(x,rho1)
grid on; yticks(0:0.2:1);
ylabel('Density Ratio')
subplot(4,1,4)
plot(x,T1)
grid on; ylim([0 1]); yticks(0:0.2:1);
ylabel('Temperature Ratio')
xlabel('X/L')
figure(2)
subplot(4,1,1)
plot(x,M2)
grid on; ylim([0 4.5]); yticks(0:1.5:4.5);
ylabel('Mach Number')
title('Flow parameters (conservative form of governing equation)')
subplot(4,1,2)
plot(x,p2)
grid on; yticks(0:0.2:1);
ylabel('Pressure Ratio')
subplot(4,1,3)
plot(x,rho2)
grid on; yticks(0:0.2:1);
ylabel('Density Ratio')
subplot(4,1,4)
plot(x,T2)
grid on; ylim([0 1]); yticks(0:0.2:1);
ylabel('Temperature Ratio')
xlabel('X/L')
figure(3)
subplot(4,1,1)
plot(1:nt,rho_t1)
hold on
plot(1:nt,rho_t2)
grid on; ylim([0 1]); yticks(0:0.2:1);
ylabel('rho/rho_0')
title('Flow parameter variation with time at the nozzle throat');
legend('Non-conservative','Conservative');
subplot(4,1,2)
plot(1:nt,T_t1)
hold on
plot(1:nt,T_t2)
grid on; yticks(0:0.2:1);
ylabel('T/T_0')
legend('Non-conservative','Conservative');
subplot(4,1,3)
plot(1:nt,p_t1)
hold on
plot(1:nt,p_t2)
grid on; yticks(0:0.2:1);
ylabel('p/p_0')
legend('Non-conservative','Conservative')
subplot(4,1,4)
plot(1:nt,M_t1)
hold on
plot(1:nt,M_t2)
grid on; yticks(0:0.2:2);
ylabel('Mach No.')
xlabel('Time Steps')
legend('Non-conservative','Conservative')
figure(4)
subplot(2,1,1)
plot(1:nt,residual1)
grid on
xlabel('Time Steps')
ylabel('Residuals')
legend('|drho/dt|_{av}','|dV/dt|_{av}','|dT/dt|_{av}')
title('Residuals: Non-conservative form')
subplot(2,1,2)
plot(1:nt,residual2)
grid on
xlabel('Time Steps')
ylabel('Residuals'); yticks(0:0.5:4);
legend('|dU_1/dt|_{av}','|dU_2/dt|_{av}','|dU_3/dt|_{av}')
title('Residuals: Conservative form')
figure(5)
hold on; grid on;
plot(x,m1(50,:),'-x')
plot(x,m1(100,:),'-d')
plot(x,m1(150,:),'-s')
plot(x,m1(200,:),'-*')
plot(x,m1(700,:),'-o')
xlabel('X/L')
ylabel('rhoAV/rho_0A_0V_0')
legend('t = 50/\t','t = 100/\t','t = 150/\t','t = 200/\t','t = 700/\t')
title('Variation of non-dimensional mass flow rate (Non-conservative form)')
figure(6)
hold on; grid on;
plot(x,m2(100,:),'-x-')
plot(x,m2(150,:),'-d-')
plot(x,m2(200,:),'-*-')
plot(x,m2(700,:),'-o-')
xlabel('X/L')
ylabel('rhoAV/rho_0A_0V_0')
legend('t = 100/\t','t = 150/\t','t = 200/\t','t = 700/\t')
title('Variation of non-dimensional mass flow rate (Conservative form)')
figure(7)
hold on; grid on;
plot(x,m1(end,:),'LineWidth',3)
plot(x,m2(end,:),'LineWidth',3)
plot(x,0.579*ones(length(x),1),'--k')
xlabel('X/L')
ylabel('rhoAV/rho_0A_0V_0')
legend('Non-conservative','Conservative','Exact')
title('Variation of non-dimensional mass flow rate (Steady State)')
2. Function for a conservative form of governing equation
% Solver for Conservative Form of Governing Equation
function [rho, V, T, p, M, m, rho_t, V_t, p_t, T_t, M_t, residual] = conservative(x, dx, A, C, nt, gamma)
% Initial conditions
rho_0 = zeros(length(x),1);
T_0 = rho_0;
V_0 = T_0;
for i=1:length(x)
if x(i)>=0 && x(i)<0.5
rho_0(i) = 1;
T_0(i) = 1;
elseif x(i)>=0.5 && x(i)<1.5
rho_0(i) = 1 - 0.366*(x(i)-0.5);
T_0(i) = 1 - 0.167*(x(i)-0.5);
elseif x(i)>=1.5 && x(i)<=3
rho_0(i) = 0.634 - 0.3879*(x(i)-1.5);
T_0(i) = 0.833 - 0.3507*(x(i)-1.5);
end
V_0(i) = 0.59/(rho_0(i)*A(i));
end
throat = find(A==1);
rho = rho_0;
T = T_0;
V = V_0;
m = zeros(nt,length(x));
U1 = rho.*A;
U2 = rho.*A.*V;
U3 = rho.*(T./(gamma-1) + 0.5*gamma*V.^2).*A;
rho_t = zeros(nt,1);
T_t = rho_t;
V_t = T_t;
p_t = V_t;
M_t = p_t;
residual = zeros(nt,3);
for k=1:nt
U1_old = U1; U2_old = U2; U3_old = U3;
dt = min(C*dx./(T.^0.5 + V));
F1 = U2;
F2 = (U2.^2)./U1+((gamma-1)/gamma)*(U3 - 0.5*gamma*(U2.^2)./U1);
F3 = gamma*U2.*U3./U1 - 0.5*gamma*(gamma-1)*(U2.^3)./(U1.^2);
% Predictor loop
for i=2:length(x)-1
dF1dx = (F1(i+1) - F1(i))/dx;
dF2dx = (F2(i+1) - F2(i))/dx;
dF3dx = (F3(i+1) - F3(i))/dx;
dAdx = (A(i+1) - A(i))/dx;
J2 = (1/gamma)*rho(i)*T(i)*dAdx;
dU1dt_p(i) = -dF1dx;
dU2dt_p(i) = -dF2dx + J2;
dU3dt_p(i) = -dF3dx;
% Intermediate update
U1(i) = U1(i) + dU1dt_p(i)*dt;
U2(i) = U2(i) + dU2dt_p(i)*dt;
U3(i) = U3(i) + dU3dt_p(i)*dt;
end
rho = U1./A;
V = U2./U1;
T = (gamma-1)*((U3./U1) - 0.5*gamma*V.^2);
F1 = U2;
F2 = (U2.^2)./U1+((gamma-1)/gamma)*(U3 - 0.5*gamma*(U2.^2)./U1);
F3 = gamma*U2.*U3./U1 - 0.5*gamma*(gamma-1)*(U2.^3)./(U1.^2);
% Corrector loop
for i=2:length(x)-1
dF1dx = (F1(i) - F1(i-1))/dx;
dF2dx = (F2(i) - F2(i-1))/dx;
dF3dx = (F3(i) - F3(i-1))/dx;
dAdx = (A(i) - A(i-1))/dx;
J2 = (1/gamma)*rho(i)*T(i)*dAdx;
dU1dt_c(i) = -dF1dx;
dU2dt_c(i) = -dF2dx + J2;
dU3dt_c(i) = -dF3dx;
end
% Average of predictor and corrector step
dU1dt = 0.5*(dU1dt_p + dU1dt_c);
dU2dt = 0.5*(dU2dt_p + dU2dt_c);
dU3dt = 0.5*(dU3dt_p + dU3dt_c);
residual(k,1) = max(abs(dU1dt));
residual(k,2) = max(abs(dU2dt));
residual(k,3) = max(abs(dU3dt));
% Final solution update
for i=2:length(x)-1
U1(i) = U1_old(i) + dU1dt(i)*dt;
U2(i) = U2_old(i) + dU2dt(i)*dt;
U3(i) = U3_old(i) + dU3dt(i)*dt;
end
% Boundary conditions
U1(1) = rho(1)*A(1);
U1(end) = 2*U1(end-1) - U1(end-2);
U2(1) = 2*U2(2) - U2(3);
U2(end) = 2*U2(end-1) - U2(end-2);
U3(1) = U1(1)*((T(1)/(gamma-1)) + 0.5*gamma*V(1)^2);
U3(end) = 2*U3(end-1) - U3(end-2);
rho = U1./A;
V = U2./U1;
T = (gamma-1)*((U3./U1) - 0.5*gamma*V.^2);
M = V./(T.^0.5);
p = rho.*T;
m(k,:) = rho.*A.*V;
rho_t(k) = rho(throat); V_t(k) = V(throat); T_t(k) = T(throat);
p_t(k) = p(throat); M_t(k) = M(throat);
end
end
3. Function for a non-conservative form of governing equation
% Solver for Non-Conservative Form of Governing Equation
function [rho, V, T, p, M, m, rho_t, V_t, p_t, T_t, M_t, residual] = nonconservative(x, dx, A, C, nt, gamma)
% Initial conditions
rho_0 = 1 - 0.3146*x;
T_0 = 1 - 0.2314*x;
V_0 = (0.1 + 1.09*x).*T_0.^0.5;
throat = find(A==1);
rho = rho_0; T = T_0; V = V_0; m = zeros(nt,length(x));
rho_t = zeros(nt,1); T_t = rho_t; V_t = T_t; p_t = V_t; M_t = p_t;
residual = zeros(nt,3);
for k=1:nt
rho_old = rho; T_old = T; V_old = V;
dt = min(C*dx./(T.^0.5 + V));
% Predictor loop
for i=2:length(x)-1
dVdx = (V(i+1) - V(i))/dx;
dlogAdx = (log(A(i+1))-log(A(i)))/dx;
drhodx = (rho(i+1) - rho(i))/dx;
dTdx = (T(i+1) - T(i))/dx;
drhodt_p(i) = -rho(i)*dVdx - rho(i)*V(i)*dlogAdx - V(i)*drhodx;
dVdt_p(i) = -V(i)*dVdx - (1/gamma)*(dTdx + (T(i)/rho(i))*drhodx);
dTdt_p(i) = -V(i)*dTdx - (gamma-1)*T(i)*(dVdx + V(i)*dlogAdx);
% Intermediate 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 loop
for i=2:length(x)-1
dVdx = (V(i) - V(i-1))/dx;
dlogAdx = (log(A(i))-log(A(i-1)))/dx;
drhodx = (rho(i) - rho(i-1))/dx;
dTdx = (T(i) - T(i-1))/dx;
drhodt_c(i) = -rho(i)*dVdx - rho(i)*V(i)*dlogAdx - V(i)*drhodx;
dVdt_c(i) = -V(i)*dVdx - (1/gamma)*(dTdx + (T(i)/rho(i))*drhodx);
dTdt_c(i) = -V(i)*dTdx - (gamma-1)*T(i)*(dVdx + V(i)*dlogAdx);
end
% Average of predictor and corrector step
drhodt = 0.5*(drhodt_p + drhodt_c);
dVdt = 0.5*(dVdt_p + dVdt_c);
dTdt = 0.5*(dTdt_p + dTdt_c);
residual(k,1) = max(abs(drhodt));
residual(k,2) = max(abs(dVdt));
residual(k,3) = max(abs(dTdt));
% Final solution update
for i=2:length(x)-1
rho(i) = rho_old(i) + drhodt(i)*dt;
V(i) = V_old(i) + dVdt(i)*dt;
T(i) = T_old(i) + dTdt(i)*dt;
end
% Boundary conditions
rho(end) = 2*rho(end-1)-rho(end-2);
V(1) = 2*V(2)-V(3); V(end) = 2*V(end-1)-V(end-2);
T(end) = 2*T(end-1)-T(end-2);
M = V./(T.^0.5);
p = rho.*T;
m(k,:) = rho.*A.*V;
rho_t(k) = rho(throat); V_t(k) = V(throat); T_t(k) = T(throat);
p_t(k) = p(throat); M_t(k) = M(throat);
end
end
Results:
Mass flow @ nx=31
mass flow @ nx=71
Conclusion:
In this challenge, we simulated the isentropic flow through a quasi-1D subsonic-supersonic nozzle. we are derived both the conservation and non-conservation forms of the governing equations and solved them using MacCormack's technique. we showed the steady-state temperature distribution for the flow-field variables and investigated the difference between the two forms of governing equations by comparing their solutions.
By using Matlab:
1. We showed the Steady-state distribution of primitive variables inside the nozzle.
2. We showed the Time-wise variation of the primitive variables.
3. We showed the Variation of Mass flow rate distribution inside the nozzle at different time steps during the time-marching process.
4. we showed the Comparison of Normalized mass flow rate distributions of both forms.