Week 7 - Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
- AIM:
- To write program in MATLAB to ‘simulate flow through a quasi 1D subsonic-supersonic nozzle’.
- Derive both the conservation and non-conservation forms of the governing equations and solve them using the MacCormack’s technique.
- Determine the steady-state temperature distribution for the flow-field variables.
- OBJECTIVE:
- Comparison of mass flow rate distribution of both conservative and nonconservative forms.
- Variation of Mass flow rate distribution inside the nozzle at different time steps during the time-marching process.
- Steady-state distribution and Time-wise variation of primitive variables.
- BODY OF THE CONTENT:
Non-Conservative Form: -
Governing Equations: -
Continuity Equation:
`(dρ^')/(dt^' )= -ρ^' (dV^')/(dx^' ) -ρ'V^' d(l n〖A^' 〗 )/dx'-V'(dρ^')/(dx^' )`
Momentum Equation:
`(dV^')/(dt^' )= -V^' (dV^')/(dx^' ) -1/γ( (dT^')/(dx^' )+T/ρ^' (dρ^')/(dx^' ) )`
Energy Equation:
`(dT^')/(dt^)`=`-V^' (dT^')/(dx^' ) -(γ-1)( (dV^')/(dx^' )+V^' d(l n〖A^' 〗 )/dx' )`
Conservative Form: -
Governing Equations: -
Continuity Equation:
`(d(ρ'V^'))/dt'= -(d(ρ'V^' A'))/dx'`
Momentum Equation:
`(d(ρ'V^' A'))/dt'+(d[ρ^' A^' V^'2+1/γ (ρ^' A^' )])/dx'= 1/γ( ρ^' (dA^')/(dx^' ))`
Energy Equation:
`(d[ρ^' (T^'/((γ-1) )+γ/2 V^'2 ) A^'] )/dt'+(d[ρ^' (T^'/((γ-1) )+γ/2 V^'2 ) (V^' A^')+ρ'A'V^'] )/dx' = 0`
- Divide the Conservative form equations for easy understanding and programming prospective is as follows
`U1= ρ'A'`
`U2= ρ'V^' A'`
`U3=[ρ^' (T^'/((γ-1) )+γ/2 V^'2 ) A^']`
`F1= U2`
`F2= [ρ^' A^' V'^2+1/γ (ρ^' A^' )`
`F3= [ρ^' (T^'/((γ-1) )+γ/2 V^'2 ) (V^' A)^'+ρ'A'V^']`
`J2= 1/γ( ρ^' (dA^')/(dx^' ))`
- Used above governing equations in program for solve the ‘quasi 1D subsonic-supersonic nozzle’ problem using MacCormack’s technique.
- We use separate functions to solve conservative and nonconservative forms.
Program Syntax:
MATLAB function for Nonconservative form
function [v,rho,t,mach,m,p] = function_nonconservative_form(n,x,dx,nt,dt_NC,a,gamma)
% calculate initial profiles
rho = 1-0.3146*x;
t = 1-0.2314*x;% t = tempreture
v = (0.1+1.09*x).*t.^0.5;
% outer time loop
for k = 1:nt
rho_old = rho;
v_old = v;
t_old = t;
% PREDECtOR method jn that we used FORWAED DjFFERENCjNG method
for j = 2:n-1
dv_dx = (v(j+1)-v(j))/dx;
drho_dx = (rho(j+1)-rho(j))/dx;
dt_dx = (t(j+1)-t(j))/dx;
dlog_a_dx = (log(a(j+1))-log(a(j)))/dx;
%continuity equation
drho_dt_precdictor(j) = -rho(j)*dv_dx-rho(j)*v(j)*dlog_a_dx-v(j)*(drho_dx);
% momemtum equatjon
dv_dt_precdictor(j) = -v(j)*dv_dx-(1/gamma)*(dt_dx + (t(j)/rho(j))*(drho_dx));
% energy equatjon
dt_dt_precdictor(j) = -v(j)*dt_dx-(gamma-1)*t(j)*(dv_dx + v(j)*(dlog_a_dx));
% solution update
% we are using this values in corrector step
rho(j) = rho(j) + drho_dt_precdictor(j)*dt_NC;
v(j) = v(j) + dv_dt_precdictor(j)*dt_NC;
t(j) = t(j) + dt_dt_precdictor(j)*dt_NC;
end
% CORRECtOR method in that we used BACKWAED DiFFERENCiNG method
for j = 2:n-1
% for continuity equation
%drho_dt_corrector=-rho(j)*(v(j)-v(j-1))/dx-rho(j)*v(j)*(log(a(j)-log(a(j-1)))/dx-v(j)*(rho(j)-rho(j-1))/dx;
% it is more lengthy so we break it is as follows and used it
drho_dx = (rho(j)-rho(j-1))/dx;
dv_dx = (v(j)-v(j-1))/dx;
dt_dx = (t(j)-t(j-1))/dx;
dlog_a_dx = (log(a(j))-log(a(j-1)))/dx;
%continuity equation
drho_dt_corrector(j) = -rho(j)*dv_dx-rho(j)*v(j)*dlog_a_dx-v(j)*drho_dx;
% momemtum equation
dv_dt_corrector(j) = -v(j)*dv_dx-(1/gamma)*((dt_dx) + (t(j)/rho(j))*drho_dx);
% energy equation
dt_dt_corrector(j) = -v(j)*dt_dx-(gamma-1)*t(j)*(dv_dx + v(j)*dlog_a_dx);
end
% compute the average tjme derjvatjve
drho_dt = 0.5*(drho_dt_precdictor + drho_dt_corrector);
dv_dt = 0.5*(dv_dt_precdictor + dv_dt_corrector);
dt_dt = 0.5*(dt_dt_precdictor + dt_dt_corrector);
% final solution update
for i = 2:n-1
rho(i) = rho_old(i) + drho_dt(i)*dt_NC;
v(i) = v_old(i) + dv_dt(i)*dt_NC;
t(i) = t_old(i) + dt_dt(i)*dt_NC;
end
% Apply the boundary condition
% inlet
v(1) = 2*v(2)-v(3);
% outlet
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);
% mASS FLOW RATE
m = rho.*a.*v;
% pressure
p = rho.*t;
% Mach Number
mach = v./sqrt(t);
figure(2)
if k == 100
plot(x,m,'b','linewidth',2)
hold on;
else if k == 200
plot(x,m,'g','linewidth',2)
else if k == 500
plot(x,m,'y','linewidth',2)
else if k == 700
plot(x,m,'r','linewidth',2)
else if k == 1000
plot(x,m,'c','linewidth',2)
grid on
hold off
legend('100dt','200dt','500dt','700dt','1000dt')
title('Mass flow rate in a nozzle wrt time step for nonconservative form');
xlabel('Non Dimensional length of the nozzle')
ylabel('mass flow rate')
end
end
end
end
end
end
end
MATLAB function for Conservative form
function[V,Rho,T,Mach,M,P] = function_conservative_form(n,x,dx,A,gamma,c,nt,L)
for i = 1:n
if x(i)>= 0 & x(i)<=L/6
Rho(i) = 1;
T(i) = 1;
end
if x(i)>=L/6 & x(i)<=L/2
Rho(i) = 1 -0.366*(x(i)-0.5);
T(i) = 1 -0.167*(x(i)-0.5);
end
if x(i)>=L/2
Rho(i) = 0.634-0.3879*(x(i)-1.5);
T(i) = 0.833-0.3507*(x(i)-1.5);
end
end
V = 0.59./(Rho.*A);
U1 = Rho.*A;
U2 = U1.*V;
U3 = U1.*((T/(gamma-1)) + (gamma*0.5*(V.^2)));
for k =1:nt
dt = min(c*(dx./(sqrt(T)+V))); % Min time step
% saving the value after each timestep
U1_old = U1;
U2_old = U2;
U3_old = U3;
F1 = U2;
F2 = ((U2.^2)./U1) + (((gamma-1)/gamma)*(U3 - ((gamma*0.5*(U2.^2))./U1)));
F3 = ((gamma*U2.*U3)./U1) - ((gamma*0.5*(gamma-1)*(U2.^3))./(U1.^2));
% Predictor Method
for j = 2:n-1
J2 = (1/gamma)*Rho(j)*T(j)*((A(j+1)-A(j))/dx);
dU1dt_p(j) = -((F1(j+1)-F1(j))/dx);
dU2dt_p(j) = -((F2(j+1)-F2(j))/dx) + J2;
dU3dt_p(j) = -((F3(j+1)-F3(j))/dx);
U1(j) = U1(j) + dU1dt_p(j)*dt;
U2(j) = U2(j) + dU2dt_p(j)*dt;
U3(j) = U3(j) + dU3dt_p(j)*dt;
% Upadated value of density, velocity and temperature after predicted values.
Rho(j) = U1(j)/A(j);
V(j) = U2(j)/U1(j);
T(j) = (gamma-1)*((U3(j)/U1(j)) - ((gamma/2)*(V(j)^2)));
end
% Update the values of F1,F2,F3 after predictor step
F1 = U2;
F2 = ((U2.^2)./U1) + (((gamma-1)/gamma)*(U3 - ((gamma*0.5*(U2.^2))./U1)));
F3 = ((gamma*U2.*U3)./U1) - ((gamma*0.5*(gamma-1)*(U2.^3))./(U1.^2));
% Corrector step
for j = 2:n-1
J2 = (1/gamma)*Rho(j)*T(j)*((A(j)-A(j-1))/dx);
dU1dt_c(j) = -((F1(j)-F1(j-1))/dx);
dU2dt_c(j) = -((F2(j)-F2(j-1))/dx) + J2;
dU3dt_c(j) = -((F3(j)-F3(j-1))/dx);
end
% Average values
dU1dt = 0.5*(dU1dt_c + dU1dt_p);
dU2dt = 0.5*(dU2dt_c + dU2dt_p);
dU3dt = 0.5*(dU3dt_c + dU3dt_p);
% Final Solution update
for i=2:n-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
% Final velocity of density, velocity and temperature.
Rho = U1./A;
V = U2./U1;
T = (gamma-1)*((U3./U1) - ((gamma/2)*(V.^2)));
% Boundary Nodes
% Inlet
U1(1) = Rho(1)*A(1);
U2(1) = 2*U2(2) - U2(3);
U3(1) = U1(1)*((T(1)/(gamma-1))+(0.5*gamma*(V(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);
M = Rho.*A.*V;
Mach = V./(sqrt(T));
P = Rho.*T;
figure(1)
ylim([0.5 0.7])
hold on
if k==100
plot(x,M,'g','linewidth',1);
else if k==200
plot(x,M,'y','linewidth',1);
else if k==500
plot(x,M,'r','linewidth',1)
else if k == 700
plot(x,M,'b','linewidth',1)
else if k == 1000
plot(x,M,'c','linewidth',1)
grid on
hold off
legend('100dt','200dt','500dt','700dt','1000dt');
title('Mass flow rate in a nozzle wrt time step for conservative form');
xlabel('Non Dimensional length of the nozzle')
ylabel('mass flow rate')
end
end
end
end
end
end
endMain Program
clear all
close all
clc
% Defining the 1D Nozzle
n = 31;
L=3;
x = linspace(0,L,n);
dx = abs(x(2)-x(1));
A = 1 +2.2*((x-(L/2)).^2);
a = A;
% Constants
gamma = 1.4; % Gas constant
c = 0.5; % Courant Number
nt = 1400; % Time steps
dt_NC = 0.0201;
% Functions
[V,Rho,T,Mach,M,P] = function_conservative_form(n,x,dx,A,gamma,c,nt,L)
[v,rho,t,mach,m,p] = function_nonconservative_form(n,x,dx,nt,dt_NC,a,gamma)
% Plots
% Mass
figure(3)
plot(x,M,'b','linewidth',2)
hold on
plot(x,m,'r','linewidth',2)
legend('Conservative','Nonconservative')
title('Mass flow rate')
xlabel('Non Dimensional length of the nozzle')
ylabel('mass flow rate')
% Pressure
figure(4)
plot(x,P,'b','linewidth',2)
hold on
plot(x,p,'r','linewidth',2)
legend('Conservative','Nonconservative')
title('Pressure')
xlabel('Non Dimensional length of the nozzle')
ylabel('pressure')
% Mach Number
figure(5)
plot(x,Mach,'b','linewidth',2)
hold on
plot(x,mach,'r','linewidth',2)
legend('Conservative','Nonconservative')
title('Mach Number')
xlabel('Non Dimensional length of the nozzle')
ylabel('Mach Number')
% Density
figure(6)
plot(x,Rho,'b','linewidth',2)
hold on
plot(x,rho,'r','linewidth',2)
legend('Conservative','Nonconservative')
title('Density')
xlabel('Non Dimensional length of the nozzle')
ylabel('density')
% Tempreture
figure(7)
plot(x,T,'b','linewidth',2)
hold on
plot(x,t,'r','linewidth',2)
legend('Conservative','Nonconservative')
title('Temperature')
xlabel('Non Dimensional length of the nozzle')
ylabel('temp')
Syntax Explanation:
- We want to solve the ‘quasi 1D subsonic-supersonic nozzle flow simulation’.
- Above we give the separate function for conservative and nonconservative form for easy understanding.
- Starting of functions give the input parameters for program i.e. Density, Temperature, Velocity.
- After that use of time loop for run the program 1400 number of time steps.
- Then use the equations above derived to solve the subsonic-supersonic nozzle flow simulation.
- Near the end of the functions plot the variation of Mass flow rate distribution inside the nozzle at different time steps during the time-marching process.
- Both the functions are used in main program to calculate the primitive variable.
- End of the main program we compare primitive variables of conservative and nonconservative in the form of graph.
- Get the outputs when run the main program is as follows
Mass flow rate in a nozzle with respect to time step for conservative form
Mass flow rate in a nozzle with respect to time step for nonconservative form
- Above graph shows that mass flow rate at different number of time steps. After 600 time steps the mass flow rate is constant.
- Following are the comparison of the conservative and nonconservative form of primitive variables.
Mass flow rate Conservative and Nonconservative Form
Pressure Conservative and Nonconservative Form
Mach Number Conservative and Nonconservative Form
Density Conservative and Nonconservative Form
Temperature Conservative and Nonconservative Form
- Above graphs of primitive variable shows that, density, pressure, temperature decreases along the length of domain and Mach number increases along the length of domain..
Grid Dependency Test:
- To conduct grid dependency test, we consider two situations with 31 grid points and 61 grid points.
- Values are shown from the grid where the throat area lies.
- Following table shows that increase the grid points increase the accuracy but it takes more time to run the simulation with more grid points.
|
Density |
Temperature |
Pressure |
Mach Number |
Mass flow rate |
|
|
Exact Analytical Solution |
0.634 |
0.833 |
0.528 |
1 |
0.579 |
|
Conservative (31 grid points) |
0.6377 |
0.8354 |
0.5327 |
0.9998 |
0.5827 |
|
Conservative (61 grid points) |
0.639 |
0.8356 |
0.5339 |
0.9979 |
0.5829 |
|
Nonconservative (31 grid points) |
0.6387 |
0.8365 |
0.5342 |
0.9994 |
0.5838 |
|
nonconservative (61 grid points) |
0.6517 |
0.8407 |
0.5479 |
0.9793 |
0.5852 |
4. Conclusion:
- Conservative form shows better mass distribution i.e. mass conservation is better than nonconservative form.
- It is very important to understand the both forms are same. Nonconservative form produces accurate result for the primitive variables and conservative form produces accurate result for the flux variables.
- Advantage of nonconservative form is on the basis of mathematical calculation effort as the conservative form requires the higher amount of work to solve the primitive variable calculation.
- Refer the grid dependency table shows that conservative form is faster than nonconservative form.