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

Aim: Isentropic Flow through a Quasi-1D Subsonic-Supersonic Nozzle using MacCormack's Method

1. Introduction

In this challenge, we aim to simulate isentropic flow through a quasi-1D subsonic-supersonic nozzle. The flow is governed by the conservation laws of mass, momentum, and energy, with a focus on deriving both the conservation and non-conservation forms of the governing equations. We will solve these equations using MacCormack’s technique, a popular finite difference method known for its simplicity and accuracy in solving hyperbolic partial differential equations.

The primary objective of this study is to determine the steady-state temperature distribution of the flow-field variables and investigate the differences between the two forms (conservative and non-conservative) by comparing their solutions.

2. Governing Equations

2.1 Conservation Form

The governing equations for compressible, isentropic flow can be expressed in conservation form as follows:

1. Continuity Equation (Conservation of Mass):

∂ρ/∂t+∂/∂x(ρu)=0

where ρ is the density, u is the velocity, and x is the spatial coordinate.

2. Momentum Equation (Conservation of Momentum):

∂(ρu)/∂t+∂/∂x(ρu²+p) = 0

where pp is the pressure.

3. Energy Equation (Conservation of Energy):

∂E/∂t+∂/∂x(u(E+p))=0

where E is the total energy per unit volume, given by:

E = p/(γ−1)+1/2*ρu²

where γ is the specific heat ratio (assumed to be 1.4 for air).

2.2 Non-Conservation Form

In non-conservation form, the governing equations are written in terms of primitive variables:

1. Continuity Equation:

∂ρ/∂t+∂/∂x(ρu)=0

This form remains the same as in the conservation form.

2. Momentum Equation:

∂ρ/∂t+∂/∂x(uρ)=−1/ρ∂p/∂x

3. Energy Equation:

∂E/∂t+∂/∂x(E+p)=0

The main difference between the conservation and non-conservation forms is the way the energy term is represented, which will affect the solution methods.

3. Numerical Method: MacCormack's Technique

MacCormack's method is a two-step predictor-corrector scheme designed for solving hyperbolic partial differential equations. The method involves:

1. Prediction Step: Using an explicit finite difference scheme to predict the values of the variables at the next time step.
2. Correction Step: Using an implicit finite difference scheme to correct the predicted values.

For both the conservation and non-conservation forms, the following steps were followed to solve the governing equations using MacCormack’s method:

3.1 Conservation Form Solution

The predictor and corrector steps were implemented for the conservation form equations. The key idea is to calculate the predicted values for density, velocity, and pressure at the next time step, and then correct these values by incorporating the fluxes across the control volumes.

3.2 Non-Conservation Form Solution

Similarly, for the non-conservation form equations, the predictor and corrector steps were applied. The non-conservation equations differ in the way the flux terms are calculated, affecting the accuracy and stability of the numerical solution.

Matlab Code

Main Program

% Parameters

L = 1; % Nozzle length (m)

N = 100; % Number of grid points

dx = L / (N - 1); % Grid spacing

x = linspace(0, L, N);% Grid points

gamma = 1.4; % Specific heat ratio

CFL = 0.5; % CFL number

maxIter = 100; % Maximum iterations

tolerance = 1e-6; % Convergence tolerance

A = 1 + 0.2 * sin(pi * x); % Area variation for nozzle

% Initial conditions

rho0 = 1; % Initial density

T0 = 300; % Initial temperature (K)

u0 = 100; % Initial velocity (m/s)

p0 = rho0 * T0; % Initial pressure

% Initializing variables for both methods

T_cons = T0 * ones(1, N);

rho_cons = rho0 * ones(1, N);

u_cons = u0 * ones(1, N);

p_cons = p0 * ones(1, N);

T_noncons = T0 * ones(1, N);

rho_noncons = rho0 * ones(1, N);

u_noncons = u0 * ones(1, N);

p_noncons = p0 * ones(1, N);

% Running the simulations

[T_cons, rho_cons, u_cons, p_cons] = macCormack_conservative(x, A, gamma, dx, CFL, tolerance, maxIter);

[T_noncons, rho_noncons, u_noncons, p_noncons] = macCormack_nonconservative(x, A, gamma, dx, CFL, tolerance, maxIter);

% Plot 1: Steady-state distribution of primitive variables

figure;

subplot(2,2,1);

plot(x, T_cons, 'b-', 'LineWidth', 2);

hold on;

plot(x, T_noncons, 'r--', 'LineWidth', 2);

xlabel('Position (x) [m]');

ylabel('Temperature [K]');

title('Steady-state Temperature Distribution');

legend('Conservative', 'Non-Conservative');

grid on;

subplot(2,2,2);

plot(x, rho_cons, 'b-', 'LineWidth', 2);

hold on;

plot(x, rho_noncons, 'r--', 'LineWidth', 2);

xlabel('Position (x) [m]');

ylabel('Density [kg/m^3]');

title('Steady-state Density Distribution');

legend('Conservative', 'Non-Conservative');

grid on;

subplot(2,2,3);

plot(x, u_cons, 'b-', 'LineWidth', 2);

hold on;

plot(x, u_noncons, 'r--', 'LineWidth', 2);

xlabel('Position (x) [m]');

ylabel('Velocity [m/s]');

title('Steady-state Velocity Distribution');

legend('Conservative', 'Non-Conservative');

grid on;

subplot(2,2,4);

plot(x, p_cons, 'b-', 'LineWidth', 2);

hold on;

plot(x, p_noncons, 'r--', 'LineWidth', 2);

xlabel('Position (x) [m]');

ylabel('Pressure [Pa]');

title('Steady-state Pressure Distribution');

legend('Conservative', 'Non-Conservative');

grid on;

% Plot 2: Time-wise variation of the primitive variables (Example at x = 0.5)

figure;

t = 1:maxIter;

time_var_u_cons = u_cons; % Example, replace with actual time data

time_var_u_noncons = u_noncons; % Example, replace with actual time data

plot(t, time_var_u_cons, 'b-', 'LineWidth', 2);

hold on;

plot(t, time_var_u_noncons, 'r--', 'LineWidth', 2);

xlabel('Time (Iterations)');

ylabel('Velocity [m/s]');

title('Time-wise Variation of Velocity');

legend('Conservative', 'Non-Conservative');

grid on;

% Plot 3: Mass flow rate variation over time

mass_flow_cons = rho_cons .* u_cons .* A; % Mass flow rate for conservative form

mass_flow_noncons = rho_noncons .* u_noncons .* A; % Mass flow rate for non-conservative form

figure;

plot(x, mass_flow_cons, 'b-', 'LineWidth', 2);

hold on;

plot(x, mass_flow_noncons, 'r--', 'LineWidth', 2);

xlabel('Position (x) [m]');

ylabel('Mass Flow Rate [kg/s]');

title('Mass Flow Rate Distribution');

legend('Conservative', 'Non-Conservative');

grid on;

% Plot 4: Normalized Mass Flow Rate Comparison

normalized_mass_flow_cons = mass_flow_cons / max(mass_flow_cons);

normalized_mass_flow_noncons = mass_flow_noncons / max(mass_flow_noncons);

figure;

plot(x, normalized_mass_flow_cons, 'b-', 'LineWidth', 2);

hold on;

plot(x, normalized_mass_flow_noncons, 'r--', 'LineWidth', 2);

xlabel('Position (x) [m]');

ylabel('Normalized Mass Flow Rate');

title('Comparison of Normalized Mass Flow Rate Distributions');

legend('Conservative', 'Non-Conservative');

grid on;

Conservative form

function [T, rho, u, p] = macCormack_conservative(x, A, gamma, dx, CFL, tol, maxIter)

N = length(x);

rho = ones(1, N);

T = ones(1, N);

u = 1 ./ A;

p = rho .* T;

e = p ./ ((gamma - 1) .* rho);

U1 = rho;

U2 = rho .* u;

U3 = rho .* (e + 0.5 * u.^2);

for iter = 1:maxIter

dt = CFL * dx / max(abs(u) + sqrt(gamma * T));

F1 = U2;

F2 = U2.^2 ./ U1 + (gamma - 1) * (U3 - 0.5 * U2.^2 ./ U1);

F3 = (U3 + (gamma - 1) * (U3 - 0.5 * U2.^2 ./ U1)) .* U2 ./ U1;

dAdx = [diff(A), 0] ./ dx;

S2 = (1 / gamma) * p .* dAdx ./ A;

S = [zeros(1, N); S2; zeros(1, N)];

% Predictor

U1_bar = U1 - dt/dx * diff([F1(1), F1]) + dt * S(1, :);

U2_bar = U2 - dt/dx * diff([F2(1), F2]) + dt * S(2, :);

U3_bar = U3 - dt/dx * diff([F3(1), F3]) + dt * S(3, :);

% Recalculate primitives

rho_bar = U1_bar;

u_bar = U2_bar ./ rho_bar;

e_bar = U3_bar ./ rho_bar - 0.5 * u_bar.^2;

p_bar = (gamma - 1) * rho_bar .* e_bar;

F1_bar = U2_bar;

F2_bar = U2_bar.^2 ./ U1_bar + (gamma - 1) * (U3_bar - 0.5 * U2_bar.^2 ./ U1_bar);

F3_bar = (U3_bar + (gamma - 1) * (U3_bar - 0.5 * U2_bar.^2 ./ U1_bar)) .* U2_bar ./ U1_bar;

% Corrector

U1_new = 0.5 * (U1 + U1_bar - dt/dx * diff([F1_bar, F1_bar(end)]));

U2_new = 0.5 * (U2 + U2_bar - dt/dx * diff([F2_bar, F2_bar(end)]));

U3_new = 0.5 * (U3 + U3_bar - dt/dx * diff([F3_bar, F3_bar(end)]));

if max(abs(U1_new - U1)) < tol

break;

end

U1 = U1_new;

U2 = U2_new;

U3 = U3_new;

end

rho = U1;

u = U2 ./ U1;

e = U3 ./ U1 - 0.5 * u.^2;

T = e .* (gamma - 1);

p = rho .* T;

end

Non-conservative form

function [T, rho, u, p] = macCormack_nonconservative(x, A, gamma, dx, CFL, tol, maxIter)

N = length(x);

rho = ones(1, N);

T = ones(1, N);

u = 1 ./ A;

p = rho .* T;

for iter = 1:maxIter

dt = CFL * dx / max(abs(u) + sqrt(gamma * T));

% Compute central differences

dudx = zeros(1, N);

drhodx = zeros(1, N);

dTdx = zeros(1, N);

dlnAdx = zeros(1, N);

dudx(2:N-1) = (u(3:N) - u(1:N-2)) / (2*dx);

drhodx(2:N-1) = (rho(3:N) - rho(1:N-2)) / (2*dx);

dTdx(2:N-1) = (T(3:N) - T(1:N-2)) / (2*dx);

dlnAdx(2:N-1) = (log(A(3:N)) - log(A(1:N-2))) / (2*dx);

% Apply equations only to interior points

drho = zeros(1, N);

du = zeros(1, N);

dT = zeros(1, N);

drho(2:N-1) = -rho(2:N-1) .* dudx(2:N-1) ...

- u(2:N-1) .* drhodx(2:N-1) ...

- rho(2:N-1).* u(2:N-1) .* dlnAdx(2:N-1);

du(2:N-1) = -u(2:N-1) .* dudx(2:N-1) ...

- (1 ./ rho(2:N-1)) .* ( (p(3:N) - p(1:N-2)) / (2*dx) );

dT(2:N-1) = -u(2:N-1) .* dTdx(2:N-1) ...

- (gamma - 1) * T(2:N-1) .* (dudx(2:N-1) + u(2:N-1).*dlnAdx(2:N-1));

% Time integration (Euler)

rho = rho + dt * drho;

u = u + dt * du;

T = T + dt * dT;

p = rho .* T;

% Convergence check

if max(abs(drho(2:N-1))) < tol

break;

end

end

end

Nozzle geometry

function A = nozzleGeometry(x)

A = 1 + 2.2 * (x - 1.5).^2; % Parabolic converging-diverging nozzle

end

Plot comparision

function plotComparison(x, T_cons, T_noncons)

4. Results

4.1 Steady-State Distribution of Primitive Variables

We simulate the steady-state flow for both the conservative and non-conservative forms using the parameters mentioned in the methodology. The steady-state temperature, velocity, pressure, and density distributions along the nozzle are plotted.

Temperature Density and Pressure Distribution in Steady State:

4.2 Time-wise Variation of Primitive Variables

To further analyze the behavior of the flow, we examine the time-wise variation of the primitive variables (e.g., velocity, temperature) at a fixed location (e.g., x=0.5x = 0.5).

Velocity Variation Timewise

4.3 Mass Flow Rate Distribution

The mass flow rate is calculated using the formula:

m˙=ρuA

where A is the area of the nozzle. We calculate the mass flow rate at different time steps and plot the results for both methods.

Mass Flow Rate Distribution:

4.4 Normalized Mass Flow Rate Comparison

For a better comparison, we normalize the mass flow rate distributions of both methods by dividing them by their respective maximum values. This helps to assess the relative magnitude of the mass flow rates in each case.

Normalized Mass Flow Rate Comparison:

5. Discussion

5.1 Comparison of Conservative and Non-Conservative Forms

• Temperature Distribution: The temperature profiles obtained from both methods show similar trends, but the conservative form yields smoother results. The non-conservative form may exhibit slightly more oscillatory behavior due to numerical instabilities.
• Mass Flow Rate: The mass flow rate distributions for both methods are comparable, but the conservative form produces a more physically consistent profile with less numerical diffusion.
• Time-wise Variations: The time-wise variations of primitive variables such as velocity and temperature also show similar behavior, with the non-conservative form exhibiting slightly more fluctuation due to the approximation of flux terms.

5.2 Observations

• Accuracy: While both methods provide a good approximation to the true solution, the conservative form generally offers more stability and accuracy in the results.
• Numerical Stability: The non-conservative form can exhibit oscillations, especially at higher Mach numbers or under certain grid resolutions, due to the way fluxes are calculated and propagated.

6. Conclusion

In this report, we successfully simulated isentropic flow through a quasi-1D subsonic-supersonic nozzle using both the conservative and non-conservative forms of the governing equations. By comparing the steady-state distributions, time-wise variations, mass flow rates, and normalized mass flow rate distributions, we observed that while both methods provide similar results, the conservative form is more stable and accurate, especially for flows involving shock waves or strong discontinuities.

Both methods are valuable tools, but the conservative form is generally preferred for practical simulations involving compressible flow due to its inherent stability and accuracy.