Week 7 - Simulation of a 1D Supersonic nozzle flow simulation using Macormack Method
Information regarding MacCormack’s technique- Maccormacks technique is a variant of the Lax-Wendorff approach but its must simpler in its application. Maccormacks technique is also explicit finite differencing technique like Lax-Wendorff approach which is second order accurate in both space and time. Maccormack technique is student friendly technique and easy for coding and understanding the program. And more importantly result obtain form maccormacks technique is satisfactory for many fluids flow application and here we used predictor step and corrected step for getting a result for both conservation and non-conservation form. In maccormack method we used forward differencing for predictor steps and backward differencing for corrector steps
A nozzle for a supersonic flow must increase in area in the flow direction, and a diffuser must decrease in area, opposite to a nozzle and diffuser for a subsonic flow. So, for a supersonic flow to develop from a reservoir where the velocity is zero, the subsonic flow must first accelerate through a converging area to a throat, followed by continued acceleration through an enlarging area.
Equations for conservation form and coding hints
In conservation approach we consider or discretise solution in three parts are solution vector, flex vector and source vector for simplification point of view
Continuity equation
`frac{partial (rhoA)}{partial t} + frac{partial (rhoAv)}{partial x} = 0`
Momentum equation
`frac{partial (rhoAv)}{partial t} + frac{partial [rhoAv^2+1/gammaPA]}{partial x} = (1/gamma)Pfrac{partialA}{partial x}`
Energy equation
`frac{partial [rho(T/(gamma-1)+gammav^2/2)A]}{partial t} + frac{partial [rhoAv(T/(gamma-1)+gammav^2/2)+rhoAv]}{partial x} = 0`
Define initial condition for `rho` and `T`. Area profile is also shown below
Solution vector
U1=`rho*a`;
U2=2*U2-U2;
U3=`U1*((T/(gamma-1))+((gamma/2)*v^2))`;
Update the old solution vector and give equation for dt
Consider u for solution vector, f for flux vector and J for source vector
Initialise flux vector for the predictor terms
Flux term
F1=U2;
F2=(U2.^2./U1)+((gamma-1)/gamma)*(U3-(gamma*U2.^2)./(2*U1));
F3=(gamma*U2.*U3./U1)-(gamma*(gamma-1)/2)*(U2.^3./U1.^2);
Source term
J1=0
`J2=(1/gamma)*rho*T*frac{dA}{dx}`;
J3=0
Update the predictor solution(differential equation) using forward difference scheme
Update the solution vector,flux vector and primitive variables using new values of predictor so as to use it in corrector step
Now, go to corrector step
Use Backward difference scheme for all variables that is used in predictor method
Then use the average of predictor and corrector step to update the derivatives
Then update the new value of solution vector at the end of each time step
Then apply boundary condition for conservative form
inlet boundary condition
U1(1)=rho(1).*a(1);
U2(1)=2*U2(2)-U2(3);
U3(1)=U1(1)*((temp(1)/(gamma-1))+((gamma/2)*V(1).^2));
Outlet boundary condition
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);
After that update the primitive variable and define pressure,Mach number and mass flow rate using the updated primitive variable
Then calculate the specific quantities at throat position.
At throat, area ratio A=1. It is at centre only.
At throat, the index is (n+1)/2;
Then plot graph for mass flow rate at different time step along the length of the nozzle
Equations for Non Conservation form and coding hints
In conservation approach we consider or discretise solution in three parts are solution vector, flex vector and source vector for simplification point of view
Continuity equation
`frac{partial rho}{partial t} = -rho frac{partial v}{partial x}-rho v frac{partial ln A}{partial x}-v frac{partial rho}{partial x}`
Momentum equation
`frac{partial v}{partial t} = -v frac{partial v}{partial x}-1/gamma (frac{partial T}{partial x}+T/rho frac{partial rho}{partial x})`
Energy equation
`frac{partial T}{partial t} = -v frac{partial T}{partial x}- (gamma-1)T (frac{partial v}{partial x}+v frac{partial lnA}{partial x})`
Define initial profile for A such that area vary as per the function of length of the nozzle
v,T,rho are primitive variable vary as per the function of length of the nozzle
Update the old values of rho,v and T at start of every time step
Use time step formula
`a=sqrt(T)`;
Use forward difference for predictor method and update the continuity,momentum and energy equation
After that update the solution variable using the predictor step
Now, go to corrector method and update theprimitive variables calculated in predictor step using backward difference scheme
Update the solution using the average of predictor and corrector method.
Update the new values of primitive variables
Define pressure,mach number and mass flow rate as a function of the updated primitive variables
Define Inlet and Outlet boundary conditions
Calculate the values of Mach number, Density, Pressure, Temperature at throat
Then plot graph for mass flow rate at different time step along the length of the nozzle
MATLAB CODE
close all; clear all;clc;
%% inputs
n=31;
x=linspace(0,3,n);
dx=x(2)-x(1);
A = 1 + 2.2*(x-1.5).^2;
%% Output
[M,m,k,rho,T,Pr,M_th,rho_th,T_th,Pr_th,residual]=macro_non_cons(x,dx,n);
[M_c,m_c,k_c,rho_c,T_c,Pr_c,M_th_c,rho_th_c,T_th_c,Pr_th_c,residual_c] = macro_cons(x,dx,n);
%% plot 3
figure(3)
subplot(4,1,1)
hold on
plot(x,M,'-');
plot(x,M_c,'-*');
legend('Non Conservative form','Conservative form');
xlabel('X'); ylabel('Mach Number');
title(sprintf('Mach Number @ Time step=%d',k));
subplot(4,1,2)
hold on
plot(x,rho,'-');
plot(x,rho_c,'-*');
legend('Non Conservative form','Conservative form');
xlabel('X'); ylabel('Density');
title(sprintf('Density @ Time step=%d',k));
subplot(4,1,3)
hold on
plot(x,T,'-');
plot(x,T_c,'-*');
legend('Non Conservative form','Conservative form');
xlabel('X'); ylabel('Temperature');
title(sprintf('Temperature @ Time step=%d',k));
subplot(4,1,4)
hold on
plot(x,Pr,'-');
plot(x,Pr_c,'-*');
legend('Non Conservative form','Conservative form');
xlabel('X'); ylabel('Pressure');
title(sprintf('Pressure @ Time step=%d',k));
%% plot 4
figure(4)
subplot(4,1,1)
hold on
plot(M_th);
plot(M_th_c);
legend('Non Conservative form','Conservative form');
xlabel('Time Step'); ylabel('Mach Number');
title(sprintf('Mach Number at centre for all Time step'));
subplot(4,1,2)
hold on
plot(rho_th);
plot(rho_th_c);
legend('Non Conservative form','Conservative form');
xlabel('Time Step'); ylabel('Density');
title(sprintf('Density at centre for all Time step'));
subplot(4,1,3)
hold on
plot(T_th);
plot(T_th_c);
legend('Non Conservative form','Conservative form');
xlabel('Time Step'); ylabel('Temperature');
title(sprintf('Temperature at centre for all Time step'));
subplot(4,1,4)
hold on
plot(Pr_th);
plot(Pr_th_c);
legend('Non Conservative form','Conservative form');
xlabel('Time Step'); ylabel('Pressure');
title(sprintf('Pressure at centre for all Time step'));
%% plot 5
figure(5)
subplot(2,1,1)
hold on
plot(residual);
title('Residual:Non Conservative form');
xlabel('Time Step'); ylabel('Residual');
legend('drho/dt','dv/dt','dT/dt');
subplot(2,1,2)
hold on
plot(residual_c);
title('Residual:Conservative form');
xlabel('Time Step'); ylabel('Residual');
legend('dU1/dt','dU2/dt','dU3/dt');
%% plot 6
figure(6)
plot(x,m);
hold on
plot(x,m_c);
xlabel('X(nozzle length)');
ylabel('Mass flow rate');
title('Mass flow rate at the end of time step in Non Conservative and Conservative form');
legend('Non Conservative form','Conservative form');Non conservative form function
function[M,m,k,rho,T,Pr,M_th,rho_th,T_th,Pr_th,residual] = macro_non_cons(x,dx,n)
%% Initial profile
rho= 1 - 0.3146*x;
T = 1 - 0.2314*x;
v = ( 0.1 + 1.09*x).*T.^(1/2);
A = 1 + 2.2*(x-1.5).^2;
gamma = 1.4;
rho_initial=rho; v_initial=v;
T_initial=T;
C=0.5;
%% Time step
nt=1400;
dt=0.01;
for k=1:nt
rho_old=rho;
v_old=v;
T_old=T;
dt=min(abs(C.*dx./(sqrt(T)+v)));
% Predictor method
for j=2:(n-1)
dvdx = (v(j+1)-v(j)) /dx;
dlogAdx = (log( A(j+1) ) - log( A(j) )) /dx;
drhodx = ( rho(j+1) - rho(j) ) /dx;
dTdx = (T(j+1)-T(j)) /dx;
% Continuity equation
drho_dt_p(j) = -rho(j)*dvdx - rho(j)*v(j)*dlogAdx - v(j)*drhodx;
% Momentum equation
dv_dt_p(j) = -v(j)*dvdx - (1/gamma)*( dTdx + (T(j)./rho(j))*drhodx );
% Energy equation
dT_dt_p(j) = -v(j)*dTdx - (gamma-1)*T(j)*(dvdx + v(j)*dlogAdx);
% Solution update
rho(j) = rho(j) + drho_dt_p(j).*dt;
v(j) = v(j) + dv_dt_p(j).*dt;
T(j) = T(j) + dT_dt_p(j).*dt;
end
% Corrector method
for j=2:(n-1)
dvdx = (v(j)-v(j-1)) /dx;
dlogAdx = (log( A(j) ) - log( A(j-1) )) /dx;
drhodx = ( rho(j) - rho(j-1) ) /dx;
dTdx = (T(j)-T(j-1)) /dx;
% Continuity equation
drho_dt_c(j) = -rho(j)*dvdx - rho(j)*v(j)*dlogAdx - v(j)*drhodx;
% Momentum equation
dv_dt_c(j) = -v(j)*dvdx - (1/gamma)*( dTdx + (T(j)./rho(j))*drhodx );
% Energy equation
dT_dt_c(j) = -v(j)*dTdx - (gamma-1)*T(j)*(dvdx + v(j)*dlogAdx);
end
% Compute average time derivative
drho_dt = 0.5*(drho_dt_p + drho_dt_c);
dv_dt = 0.5*(dv_dt_p + dv_dt_c);
dT_dt = 0.5*(dT_dt_p + dT_dt_c);
residual(k,1) = max(abs(drho_dt));
residual(k,2) = max(abs(dv_dt));
residual(k,3) = max(abs(dT_dt));
%% Update for the new values of primitive variable
for i=2:(n-1)
rho(i) = rho_old(i) + drho_dt(i)*dt;
v(i) = v_old(i) + dv_dt(i)*dt;
T(i) = T_old(i) + dT_dt(i)*dt;
end
%% Define pressure,mach number and mass flow rate
Pr = rho.*T;
M = v./sqrt(T);
m = rho.*A.*v;
%% Inlet BC
v(1) = 2*v(2) - v(3);
%% Outlet BC
rho(n) = 2*rho(n-1) - rho(n-2);
v(n) = 2*v(n-1) - v(n-2);
T(n) = 2*T(n-1) - T(n-2);
%% Calculating the values of quantities at throat
M_th(k)=M((n+1)/2);
rho_th(k)=rho((n+1)/2);
Pr_th(k)=Pr((n+1)/2);
T_th(k)=T((n+1)/2);
%% comparison of mass flow rates at different timesteps
figure(1)
grid on
if k==50
plot(x,m,'g')
end
hold on
if k== 100
plot(x,m,'r')
end
if k == 150
plot(x,m,'b')
end
if k==200
plot(x,m,'k')
end
if k==700
plot(x,m,'m')
end
end
title('Mass flow rates at different time steps- Non Conservative form');
legend('50','100','150','200','700');
xlabel('X(nozzle length)'); ylabel('Mass flow rate');
endConservative form function
function[M,m,k,rho,T,Pr,M_th,rho_th,T_th,Pr_th,residual] = macro_cons(x,dx,n)
%% Time step
nt=1400;
dt=0.01;
%% Initial condition for rho and T
for i=1:length(x)
if x(i)>=0 && x(i)<=5
rho(i) = 1;
T(i) = 1;
end
if x(i)>= 0.5 && x(i)<= 1.5
rho(i)=1-0.366*(x(i)-0.5);
T(i)=1-0.167*(x(i)-0.5);
end
if x(i)>=1.5 && x(i)<= 3.5
rho(i)=0.634-0.3879*(x(i)-1.5);
T(i)=0.833-0.3507*(x(i)-1.5);
end
end
%% Initial profile
A = 1 + 2.2*(x-1.5).^2;
gamma = 1.4;
%% Initializing velocity and pressure
v=0.59./(rho.*A);
Pr=rho.*T;
rho_initial=rho; v_initial=v;
T_initial=T;
C=0.5;
%% Initializing the solution vector
U1=rho.*A;
U2=rho.*A.*v;
U3=rho.*A.*( (T./(gamma-1))+(gamma/2)*v.^2);
for k=1:nt
%% updating old solution vectors
U1_old=U1;
U2_old=U2;
U3_old=U3;
%% time steps for conservative form
dt=min(abs((C*dx)./(sqrt(T)+v)));
%% predictor
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);
for j=2:(n-1)
%predictor step
J2(j)=(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(j);
dU3dt_p(j)=-(F3(j+1)-F3(j))/dx;
%updating the new values of solution vectors
U1(j)=U1(j)+dU1dt_p(j)*dt;
U2(j)=U2(j)+dU2dt_p(j)*dt;
U3(j)=U3(j)+dU3dt_p(j)*dt;
end
%% update the primitive variables and flux vectors
rho=U1./A;
v=U2./U1;
T=(gamma-1)*(U3./U1- (gamma/2)*v.^2);
Pr=rho.*T;
F1=U2;
F2=(U2.^2./U1)+((gamma-1)/gamma)*(U3-(gamma*U2.^2)./(2*U1));
F3=(gamma*U2.*U3./U1)-(gamma*(gamma-1)/2)*(U2.^3./U1.^2);
%%
for j=2:(n-1)
%% corrector step
J2(j)=(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(j);
dU3dt_c(j)=-(F3(j)-F3(j-1))/dx;
end
%% Updating the average value
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));
%% Update new value
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
%% applying BC
%% inlet BC
U1(1)=rho(1).*A(1);
U2(1)=2*U2(2)-U2(3);
U3(1)=U1(1).*((T(1)/(gamma-1))+((gamma/2)*v(1).^2));
%% Outlet BC
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);
%% updation of the primitive variables
rho = U1./A;
v = U2./U1;
T=(gamma-1)*(U3./U1-(gamma/2)*v.^2);
%% defining values for Mach number ,flow rate and pressure
M = v./sqrt(T);
m =rho.*A.*v;
Pr=rho.*T;
%% calculating the values of specific quantities at throat position
M_th(k)=M((n+1)/2);
rho_th(k)=rho((n+1)/2);
Pr_th(k)=Pr((n+1)/2);
T_th(k)=T((n+1)/2);
%% comparison mass flow rates at different timesteps
figure(2);
hold on;
grid on;
if k== 1
plot(x,m,'b')
end
if k== 100
plot(x,m,'g')
end
if k== 200
plot(x,m,'k')
end
if k== 700
plot(x,m,'m')
end
end
title('Mass flow rates at different time steps- Conservative form');
legend('1','100','200','700');
xlabel('X(nozzle length)'); ylabel('mass flow rate');
endOUTPUT
1) At n=31
a)Mass flow rates at different time steps- Non Conservative form
b)Mass flow rates at different time steps- Conservative form
c) Mach number, Density, Pressure, Temperature variation along X axis after the end of time step
d) Mach number, Density, Pressure, Temperature at throat after each timestep
e) Residuals at exit at every time step
f) Mass flow rates at different time steps- Non Conservative form
CONCLUSION
From figure a and figure b, it is noticed that the mass flow rate attains steady state nearly at 700 time step. The variation gets reduced from 0 to 700 time step.
From figure c, it is observed that the values of Mach Number, Pressure, Density and Temperature are almost same from both conservative and non conservative form of governing equation.
From figure d, it is observed that the values of Mach Number, Pressure, Density and Temperature attained steady state quicker in conservative form than that of in non conservative form. The variation is small
Figure e is just a residual plot at the nozzle exit at every time step
Figure f shows that the mass flow rate at the end of simulation is nearly constant in conservative form and it is varying in between 0.585 to 0.595 in non - conservative form
2) At n=91
a)Mass flow rates at different time steps- Non Conservative form
Use 1400 instead of 700
b)Mass flow rates at different time steps- Conservative form
Use 1400 instead of 700
c) Mach number, Density, Pressure, Temperature variation along X axis after the end of time step
d) Mach number, Density, Pressure, Temperature at throat after each timestep
e) Residuals at exit at every time step
f) Mass flow rates at different time steps- Non Conservative form
CONCLUSION
From figure a and figure b, it is found that as the mesh size reduced, it takes more time to attain steady state. At 700 time step, the solution has not reached steady state, when compared to case at n=31.
From figure d, it is found that, it takes more time for non conservative form to attain steady state at throttle area, but for conservative form, it attains steady state faster irrespective of mesh size
From figure f, it is found that for non conservative form, there is not much variation in mass flow rate at outlet when compared to case at n=31. For conservative form, the value is not changing with the length of the nozzle