Week 7 - Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
Simulation of 1D supersonic nozzel flow using Macormack method
Objectives:
1. Steady-state distribution of premitive variable inside the nozzel.
2. Time-wise distribution of premitive variables.
3. Variation of mass flow rate distribution inside the nozzel at different time steps during the time marching process.
4. Comparison of normalised mass flow rate distribution of both forms.
We will be simulating the isentropic flow through a Quasi 1D subsonic-supersonic nozzel. For this simulation we will be using both forms of governing equations, namely non-conservative and conservative forms. We will first go through the idea behid each form and then we will see the code to simulate the quasi 1D nozzel flow.
The nozzel model has a convergent section as the inlet and the divergent section as the outlet, the inlet is attached to a resiorvior so that we can assume that teh velocity of flow at the inlet is approaching zero and the area tends to infinity the pressure and temperature at the inlet are the total pressure and total temperature for the resiorvior. The mach number for the inlet is less than one that is the flow is subsonic. The outlet of the nozzel is supersonic so the mach number will be greater than one. Also at the throat of the nozzel we assume that the flow has a mach number value equal to unity.
Consider the following image for the quasi 1D supersonic nozzel flow
Non-Conservative Form:
Now we will see the non conservative form of governing equation that are associated with this flow.
Continuity equation:
`(del rho)/(del t) = - rho (del v)/(del x) - rho v (del ln a)/(delx) - v(del rho)/(del x)`
Momentum equation:
`(del v)/(del t) = - v(del v)/(del x) - 1/gamma( (del T)/(del x) + tT/rho (del rho)/ (del x))`
Enegry equation:
`(del T )/ (del t) = - v (del T) /(del x) - (gamma -1) T ((delv)/(del x) + v(del (ln a))/(delx))`
here all the variable are in nondimensionalized forms.
We willbe using macormack method for our numerical simulations
the time step is taken as: `dt = min( CFL * dx/((sqrtT) + V))`
Now we will see the function code for the non conservative form:
function [rho_n, v_n, t_n, p_n, M_n, mfr_n, rho_th_n, v_th_n, t_th_n, p_th_n, M_th_n, mfr_th_n] = non_conservative(n, nt,x, gamma, CFL, a, dx, throat)
%% function for Quasi 1D nozzel flow using non-conservative form
% calculate initial profiles for non conservative form
rho_n = 1 - 0.3146*x; % density
t_n = 1- 0.2314*x; % temperature
v_n = (0.1 + 1.09*x).*t_n.^0.5; % velocity
% time loop
for k = 1:nt
% minimum time step value among all the nodes approach should be
% employed.
dt = min(CFL*dx./((t_n.^0.5)+v_n));
rho_old = rho_n;
v_old = v_n;
t_old = t_n;
%% predictor method
for j = 2:n-1
% spatial derivativs
dvdx = (v_n(j+1) - v_n(j))/dx;
drhodx = (rho_n(j+1) - rho_n(j))/dx;
dlogadx = (log(a(j+1)) - log(a(j)))/dx;
dtdx = (t_n(j+1) - t_n(j))/dx;
% predictive step
% continuity equation
drho_dt_p(j) = - rho_n(j)*dvdx - rho_n(j)*v_n(j)*dlogadx - v_n(j)*drhodx;
% momentum equation
dv_dt_p(j) = - v_n(j)*dvdx - (1/gamma)*(dtdx + (t_n(j)/rho_n(j)*drhodx));
% energy equation
dt_dt_p(j) = -v_n(j)*dtdx - (gamma -1)*t_n(j)*(dvdx + v_n(j)*dlogadx);
% solution update
v_n(j) = v_n(j) + dv_dt_p(j) * dt;
rho_n(j) = rho_n(j) + drho_dt_p(j) * dt;
t_n(j) = t_n(j) + dt_dt_p(j) * dt;
end
%% corrector method
for j = 2:n-1
% spatial derivatives
dvdx = (v_n(j) - v_n(j-1))/dx;
drhodx = (rho_n(j) - rho_n(j-1))/dx;
dlogadx = (log(a(j)) - log(a(j-1)))/dx;
dtdx = (t_n(j) - t_n(j-1))/dx;
% corrector step
% continuity equation
drho_dt_c(j) = - rho_n(j)*dvdx - rho_n(j)*v_n(j)*dlogadx - v_n(j)*drhodx;
% momentum equation
dv_dt_c(j) = - v_n(j)*dvdx - (1/gamma)*(dtdx + (t_n(j)/rho_n(j)*drhodx));
% energy equation
dt_dt_c(j) = -v_n(j)*dtdx - (gamma -1)*t_n(j)*(dvdx + v_n(j)*dlogadx);
end
%% compute the average time derivative
dv_dt = 0.5*(dv_dt_p + dv_dt_c);
drho_dt = 0.5*(drho_dt_p + drho_dt_c);
dt_dt = 0.5*(dt_dt_p + dt_dt_c);
%% final solution update
for i = 2:n-1
v_n(i) = v_old(i) + dv_dt(i) * dt;
rho_n(i) = rho_old(i) + drho_dt(i) * dt;
t_n(i) = t_old(i) + dt_dt(i) * dt;
end
%% apply boundary conditions
% inlet
v_n(1) = 2* v_n(2) - v_n(3);
% outlet
v_n(n) = 2*v_n(n-1) - v_n(n-2);
rho_n(n) = 2*rho_n(n-1) - rho_n(n-2);
t_n(n) = 2*t_n(n-1) - t_n(n-2);
% other properties
p_n = rho_n.*t_n; % pressure
M_n = v_n./(t_n.^0.5); % Mach number
mfr_n=rho_n.*v_n.*a; % mass flow rate
% calculating the values for throat
rho_th_n(k) = rho_n(throat);
v_th_n(k) = v_n(throat);
t_th_n(k) = t_n(throat);
p_th_n(k) = p_n(throat);
M_th_n(k) = M_n(throat);
mfr_th_n(k) = mfr_n(throat);
%% plotting the mass flow rate for various time steps
figure(1)
hold on
title(' distribution of mass flow rate over the length of nozzel at various time steps')
if k == 1 % at time step = 1
plot(x,mfr_n,'--','linewidth',2)
elseif k == 10 % at time step = 10
plot(x,mfr_n,'b','linewidth',2)
elseif k == 50 % st time step = 50
plot(x,mfr_n,'r','linewidth',2)
elseif k == 100 % at time step = 100
plot(x,mfr_n,'y','linewidth',2)
elseif k == 200 % at time step = 200
plot(x,mfr_n,'c','linewidth',2)
elseif k == 700 % at time step = 700
plot(x,mfr_n,'g','linewidth',2)
legend('1','10','50','100','200','700')
xlabel('nondimensional distance through nozzel')
ylabel('mass flow rate')
end
end
%% plottting the timewise distribution of the variables at the throat
figure(2)
subplot(5,1,1)
plot(rho_th_n)
ylabel('density')
title('timewise distribution of premitive variables at the throat')
subplot(5,1,2)
plot(v_th_n)
ylabel('velocity')
subplot(5,1,3)
plot(t_th_n)
ylabel('temperature')
subplot(5,1,4)
plot(p_th_n)
ylabel('pressure')
subplot(5,1,5)
plot(M_th_n)
ylabel('mach number')
xlabel('number of time steps')
%% ploting the distribution of variables over the length of the nozzel
figure(3)
subplot(5,1,1)
plot(x,v_n)
ylabel('velocity')
title(' distribution of premitive variables over the length of nozzel')
hold on
subplot(5,1,2)
plot(x,t_n)
ylabel('temperature')
hold on
subplot(5,1,3)
plot(x,rho_n)
ylabel('density')
hold on
subplot(5,1,4)
plot(x,p_n)
ylabel('pressure')
hold on
subplot(5,1,5)
plot(x,M_n)
ylabel('Mach number')
xlabel('X')
endConstervative Form:
We will now see the conservative form of governing equations associated with the flow
Continuity equation:
`(del rho A)/ (delt) + (del rho A v )/ (delx) = 0`
Momentum equation:
`(del rho A v)/ (del t) + (del rho A v^2 + 1/gamma p A)/(del x) = 1/ gamma p (del A)/ (delx)`
Energy equation:
`(del (rho(e/ (gamma -1) +gamma/2 v^2)A))/ (del t) + (del (rho(e/(gamma-1) +gamma/2 v^2) A v + p A v))/(delx) = 0`
all the above equations are in the form of non dimensioinalised variables.
We will be converting our above equation in the form of solution vector , flux vectors and source terms
solution vectors:
`U1 = rho A`
`U2 = rho A v`
`U3 = rho ( e/ (gamma -1) + gamma/2 v^2) A`
flux vectors:
`F1 = rho A v`
`F2 = rho A v^2 + 1/ gamma p A`
`F3 = rho ( e / (gamma-1) +gamma/2 v^2) v A + p A v`
source terms:
`J2 = 1/ gamma p (del A)/ (del x)`
now the after substituting the values for the solution vectors , flux vectors and source terms our governing equation will look like the following
Continuity equation:
`(del U1)/ (del t) = - (del F1)/ (del x)`
Momentum equation:
`(del U2)/(delt) = - (delF2)/ (delx) + J2`
Energy equation:
`(del U3)/(delt) = -(delF3)/ (delx)`
we will use this above form of governing equation for our macormack method in our code
Now we will see the function code for our conservative form:
function [ rho_c, v_c, t_c, p_c, M_c, mfr_c, rho_th_c, v_th_c, t_th_c, p_th_c, M_th_c, mfr_th_c] = conservative(n, nt,x, gamma, CFL , a , dx, throat)
%% Function for the quasi 1D nozzel flow using conservative form
% calculate initial profile for conservative form
n_x_1 = ceil(0.5*(n-1)/3.5 + 1); % number of node at x = 0.5
n_x_2 = ceil(1.5*(n-1)/3.5 + 1); % number of node at x = 1.5
% density
rho_c(1:n_x_1) = 1.0;
rho_c(n_x_1:n_x_2) = 1.0 - 0.366*(x(n_x_1:n_x_2)-0.5);
rho_c(n_x_2:n) = 0.634 - 0.3879*(x(n_x_2:n)-1.5);
% temperature
t_c(1:n_x_1) = 1.0;
t_c(n_x_1:n_x_2) = 1.0 - 0.167*(x(n_x_1:n_x_2)-0.5);
t_c(n_x_2:n) = 0.833 - 0.3507*(x(n_x_2:n)-1.5);
% velocity
v_c = 0.59./(rho_c.*a);
% converting the premitive variables into solution variables
% solution vectors
U1 = rho_c.*a; % continuity
U2 = rho_c.*a.*v_c; % momentum
U3 = rho_c.*(t_c./(gamma-1) + gamma/2*v_c.^2).*a; % energy
% time loop
for k = 1:nt
% minimum time step value among all the nodes approach should be
% employed.
dt = min(CFL*dx./(sqrt(t_c(2:n-1)) + v_c(2:n-1)));
% saving the solution vectors
U1_old = U1;
U2_old = U2;
U3_old = U3;
% flux vectors
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;
%% perdictor method
for i = 2:n-1
% source term
%J2 = (gamma -1)/gamma * (U3 - gamma/2 * U2.^2./U1) * dadx;
J2_p(i) = (1/gamma)*rho_c(i)*t_c(i)*((a(i+1) - a(i))/dx);
% predictor step
% forward diferencing
dU1dt_p(i) = - (F1(i+1) - F1(i))/dx; % continuity
dU2dt_p(i) = - (F2(i+1) - F2(i))/dx + J2_p(i); % momentum
dU3dt_p(i) = - (F3(i+1) - F3(i))/dx; % energy
% updating the solution variables
U1(i) = U1_old(i) + (dU1dt_p(i)*dt);
U2(i) = U2_old(i) + (dU2dt_p(i)*dt);
U3(i) = U3_old(i) + (dU3dt_p(i)*dt);
end
%% decoding the primitive variables
rho_c_p = U1./a;
v_c_p = U2./U1;
t_c_p = (gamma -1)*(U3./U1 - (gamma/2) * ((v_c_p).^2));
F1_c = U2;
F2_c = U2.^2./U1 + (gamma -1)/gamma * (U3 - gamma/2 * U2.^2./U1);
F3_c = gamma*U2.*U3./U1 - gamma*(gamma-1)/2* U2.^3./U1.^2;
%% corrector method
for i = 2:n-1
% source terms
J2_c(i) = (1/gamma)*rho_c_p(i)*t_c_p(i)*((a(i)-a(i-1))/dx);
% corrector step
% backward differencing
dU1dt_c(i) = - (F1_c(i) - F1_c(i-1))/dx; % continutiy
dU2dt_c(i) = - (F2_c(i) - F2_c(i-1))/dx + J2_c(i); % momentum
dU3dt_c(i) = - (F3_c(i) - F3_c(i-1))/dx; % energy
end
% average time derivatives
dU1dt_avg = 0.5*(dU1dt_p+ dU1dt_c);
dU2dt_avg = 0.5*(dU2dt_p+ dU2dt_c);
dU3dt_avg = 0.5*(dU3dt_p+ dU3dt_c);
for i = 2:n-1
% final updated values
U1(i) = U1_old(i) + dU1dt_avg(i)*dt; % continuity
U2(i) = U2_old(i) + dU2dt_avg(i)*dt; % momentum
U3(i) = U3_old(i) + dU3dt_avg(i)*dt; % energy
end
%% apply boundary conditions
% inlet
U1(1)=rho_c_p(1).*a(1);
U2(1) =2*U2(2) - U2(3);
v_c_p(1) = U2(1)/U1(1);
U3(1)= U1(1)*(t_c_p(1)/(gamma-1)+(gamma/2)*v_c_p(1).^2);
% outlet
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 updated primitive variables
rho_c = U1./a ;
v_c = U2./U1;
t_c = (gamma - 1)* (U3/U1 - (gamma/2)*v_c.^2);
% other premitive variables
p_c = rho_c.*t_c; % pressure
M_c = v_c./(t_c.^0.5); % Mach number
mfr_c =rho_c.*v_c.*a; % mass flow rate
% values at throat
rho_th_c(k) = rho_c(throat);
v_th_c(k) = v_c(throat);
t_th_c(k) = t_c(throat);
M_th_c(k) = M_c(throat);
p_th_c(k) = p_c(throat);
mfr_th_c(k) = mfr_c(throat);
%% plotting the mass flow rate for various time steps
figure(4)
hold on
title(' distribution of mass flow rate over the length of nozzel at various time steps')
if k == 1
plot(x,mfr_c,'--','linewidth',2)
elseif k == 10
plot(x,mfr_c,'b','linewidth',2)
elseif k == 50
plot(x,mfr_c,'r','linewidth',2)
elseif k == 100
plot(x,mfr_c,'y','linewidth',2)
elseif k == 200
plot(x,mfr_c,'c','linewidth',2)
elseif k == 700
plot(x,mfr_c,'g','linewidth',2)
legend('1','10','50','100','200','700')
xlabel('nondimensional distance through nozzel')
ylabel('mass flow rate')
end
end
%% plottting the timewise distribution of the variables at the throat
figure(5)
subplot(5,1,1)
plot(rho_th_c)
ylabel('density')
title('timewise distribution of primitive variables at the throat')
subplot(5,1,2)
plot(v_th_c)
ylabel('velocity')
subplot(5,1,3)
plot(t_th_c)
ylabel('temperature')
subplot(5,1,4)
plot(p_th_c)
ylabel('pressure')
subplot(5,1,5)
plot(M_th_c)
ylabel('mach number')
xlabel('number of time steps')
%% ploting the distribution of variables over the length of the nozzel
figure(6)
subplot(5,1,1)
plot(x,v_c)
ylabel('velocity')
title(' distribution of premitive variables over the length of nozzel')
hold on
subplot(5,1,2)
plot(x,rho_c)
ylabel('density')
hold on
subplot(5,1,3)
plot(x,t_c)
ylabel('temperature')
hold on
subplot(5,1,4)
plot(x,p_c)
ylabel('pressure')
hold on
subplot(5,1,5)
plot(x,M_c)
ylabel('mach number')
xlabel('X')
endNow using both the functions for Non-consesrvative and conservative forms in our main program we will see the results for a 1D nozzel geometry with 31 nodes.
clear all
close all
clc
%inputs
n = 31; % number of nodes
x = linspace(0,3,n);
dx = x(2)-x(1); % element length
gamma = 1.4;
% area
a = 1 + 2.2*(x - 1.5).^2;
%CFL criteria
CFL = 0.5;
% throat of the nozzel
throat = find(a==1);
% number of time steps
nt = 1400;
%% using the functions for non-conservative and conservative form of governing equations
% for non-conservative form
[rho_n, v_n, t_n, p_n, M_n, mfr_n, rho_th_n, v_th_n, t_th_n, p_th_n, M_th_n, mfr_th_n] = non_conservative(n,nt,x,gamma,CFL,a,dx,throat);
% for conservative form
[rho_c, v_c, t_c, p_c, M_c, mfr_c, rho_th_c, v_th_c, t_th_c, p_th_c, M_th_c, mfr_th_c] = conservative(n,nt,x,gamma,CFL,a,dx,throat);
%% plotting the mass flow rate comparison at the steady state for both forms
figure(7)
plot(x,mfr_c,'linewidth',2)
hold on
plot(x,mfr_n,'linewidth',2)
legend('conservative','non-conservative')
xlabel('non dimensional length of nozzel')
ylabel('non dimensional mass flow rate')
title('comparison of normalised mass-flow-rate distributions for both forms')
OUTPUTS:
The results for the main program is given below:
Mass flow rate at various time steps for non conservative form:
Time wise distribution of variables at the throat for non conservative form:
here we can see athat after 800 timesteps the values of the variables does not change.
Distribution of properties across the length of the nozzel for non conservative form:
Mass flow rate at various time steps for conservative form:
Time wise distribution of variables at the throat for conservative form:
here we can see athat after 800 timesteps the values of the variables does not change.
Distribution of properties across the length of the nozzel for conservative form:
Also we will see the comparison of mass flow rate for both forms
Comparison of mass flow rates:
From the above plot we can see that the mass flow rate for the conservative form remains almost constant throughout the length of the nozzzel where as in the non conservative form the mass flow rate fluctuates more. in the Non-conservative forms the residuals and the round off error are very less as compared to the conservative forms. This is due to the conversion from the primitive variables to solution vectors, and back then to primitive variables in the conservative form.
Grid Independence test:
Now we will perform the grid independence test for both forms:
we will now use 61 nodes in our simulation
the code in the main program will have 61 as the number of nodes and all the other things will remian same for the code.
similarly ,we will also use 91 grid point to perform the test.
the code for grid indepndence test is given below:
clear all
close all
clc
%inputs
for n = 31; % number of nodes
x = linspace(0,3,n);
dx = x(2)-x(1); % element length
gamma = 1.4;
% area
a = 1 + 2.2*(x - 1.5).^2;
%CFL criteria
CFL = 0.5;
% throat of the nozzel
throat = find(a==1);
% number of time steps
nt = 1400;
%% using the functions for non-conservative and conservative form of governing equations
% for non-conservative form
[rho_n, v_n, t_n, p_n, M_n, mfr_n, rho_th_n, v_th_n, t_th_n, p_th_n, M_th_n, mfr_th_n] = non_conservative(n,nt,x,gamma,CFL,a,dx,throat);
% for conservative form
[rho_c, v_c, t_c, p_c, M_c, mfr_c, rho_th_c, v_th_c, t_th_c, p_th_c, M_th_c, mfr_th_c] = conservative(n,nt,x,gamma,CFL,a,dx,throat);
rho_th_n_gr(1) = rho_th_n(1400)
v_th_n_gr(1) = v_th_n(1400)
t_th_n_gr(1) = t_th_n(1400)
rho_th_c_gr(1) = rho_th_c(1400)
v_th_c_gr(1) = v_th_c(1400)
t_th_c_gr(1) = t_th_c(1400)
end
%inputs
for n = 61; % number of nodes
x = linspace(0,3,n);
dx = x(2)-x(1); % element length
gamma = 1.4;
% area
a = 1 + 2.2*(x - 1.5).^2;
%CFL criteria
CFL = 0.5;
% throat of the nozzel
throat = find(a==1);
% number of time steps
nt = 1400;
%% using the functions for non-conservative and conservative form of governing equations
% for non-conservative form
[rho_n, v_n, t_n, p_n, M_n, mfr_n, rho_th_n, v_th_n, t_th_n, p_th_n, M_th_n, mfr_th_n] = non_conservative(n,nt,x,gamma,CFL,a,dx,throat);
% for conservative form
[rho_c, v_c, t_c, p_c, M_c, mfr_c, rho_th_c, v_th_c, t_th_c, p_th_c, M_th_c, mfr_th_c] = conservative(n,nt,x,gamma,CFL,a,dx,throat);
rho_th_n_gr(2) = rho_th_n(1400)
v_th_n_gr(2) = v_th_n(1400)
t_th_n_gr(2) = t_th_n(1400)
rho_th_c_gr(2) = rho_th_c(1400)
v_th_c_gr(2) = v_th_c(1400)
t_th_c_gr(2) = t_th_c(1400)
end
%inputs
for n = 91; % number of nodes
x = linspace(0,3,n);
dx = x(2)-x(1); % element length
gamma = 1.4;
% area
a = 1 + 2.2*(x - 1.5).^2;
%CFL criteria
CFL = 0.5;
% throat of the nozzel
throat = find(a==1);
% number of time steps
nt = 1400;
%% using the functions for non-conservative and conservative form of governing equations
% for non-conservative form
[rho_n, v_n, t_n, p_n, M_n, mfr_n, rho_th_n, v_th_n, t_th_n, p_th_n, M_th_n, mfr_th_n] = non_conservative(n,nt,x,gamma,CFL,a,dx,throat);
% for conservative form
[rho_c, v_c, t_c, p_c, M_c, mfr_c, rho_th_c, v_th_c, t_th_c, p_th_c, M_th_c, mfr_th_c] = conservative(n,nt,x,gamma,CFL,a,dx,throat);
rho_th_n_gr(3) = rho_th_n(1400)
v_th_n_gr(3) = v_th_n(1400)
t_th_n_gr(3) = t_th_n(1400)
rho_th_c_gr(3) = rho_th_c(1400)
v_th_c_gr(3) = v_th_c(1400)
t_th_c_gr(3) = t_th_c(1400)
end
%% plotting the mass flow rate comparison at the steady state for both forms
figure(7)
plot(x,mfr_c,'linewidth',2)
hold on
plot(x,mfr_n,'linewidth',2)
legend('conservative','non-conservative')
xlabel('non dimensional length of nozzel')
ylabel('non dimensional mass flow rate')
title('comparison of normalised mass-flow-rate distributions for both forms')
From the simulation we see that the result
we will first compare the results of the non conservative from
now the results of the conservative form:
from the results we can see that there is not much change in the values of the variables so our numerical simulation is grid independent and the grid is good for teh results.
Also we can see that there is less change in the values of the non conservational form as compared to the comservational form.