Week 7 - Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
Aim: Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method.
Objectives:
1.Simulate the isentropic flow through a quasi 1D subsonic-supersonic nozzle.
2.Derive both the conservation and non-conservation forms of the governing equations and sovle them using the MacCormack's technique.
3. 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.
Description: We consider steady,isentropic flow through convergent nozzle. The flow to the inlet to the nozzle comes from the nozzle where reservoir pressure temperature and density are `p_0,T_0 and rho_0`respectilvely. Area of reservoir is considered to be very large as a result the velocity is considered as zero therefore `p_0,T_0 and rho_0` becomes satgnation pressure temperature and density. The flow enters the nozzle at subsonis smeed i.e M < 1 and it reaches supersoniv values in the divergent section of the nozzle, At the throat where Area is 1 we reach sonic value that is M = 1.
Governing equation describing the flow are:
These cannot be directly applied to mackormack methof for CFD simulation, therefore first we will derive the equations suitable for time marching technique fro mackormack method.
Setup: we derive the governing equations by considering small control volume region inside the nozzle.
By applying the physical princeples of the flow we obtain conservative and non conservative form of governing equations, We obtain non conservative form of governing equation when the control volume is fixed in sapce, If the control volume is moving we obtain conservative form of governing equation.
The governing equations for quasi one dimensional steady and isentropic flow are:
Non conservative forms of governing equations:
Continuity equation:
`(del(rho^'))/(delt^')=-rho^'(delV^')/(delx^')-rho^'V^'(del(lnA^'))/(delx^')-V^'(delrho^')/(delx^')`
Momentum equation:
`(delV^')/(delt^') = -V^' (delV^')/(delx^') - 1/gamma((delT^')/(delx^')+T^'/rho^'(delrho^')/(delx^'))`
Energy equation:
`(delT^')/(delt^') = -V^'(delT^')/(delx^') -(gamma-1)T^'[(delV^')/(delx^') +V^'(del(lnA^'))/(delx^') ]`
Conservative form of governing equation:
Continuity equation:
`(del(rho^'A^'))/(delt^')+(del(rho^'A^'V^'))/(delx^') = 0`
Momentum equation:
`(del(rho^'A^'V^'))/(delt^') + (del[rho^'A^'(V^')^2+(1/gamma)rho^'A^'])/(delx^')= (1/gamma)rho^'((delA^')/(delx^'))`
Energy equation:
`(del[rho^'(e^'/(gamma-1)+(gamma/2)(V^')^2)A^'])/(delt^') + (del[rho^'(e^'/(gamma-1)+(gamma/2)(V^')^2)A^'V^'+rho^'A^'V^'])/(delx^')=0`
where `gamma = 1.4`
We apply macormack method to these equations.
Macormack method: Macormack method is also an explicit finite difference technique which is second order accurate in both saoce and time.
let us explain macormack method:
value of f at next time step by macormack method is given by:
`f_(i,j)^(t+Deltat) = f_(i,j)^t + ((delf)/(delt))_(avg) Deltat` ---------- equation 1
`((delf)/(delt))_(avg)`is a representative mean value between t and `t+Deltat`.
macormack method is done in two steps:
1. Predictor step.
2. Corrector step.
Predictor step: In predictorstep we find the value of `((delf)/(delt))`at time t by forward difference we than update the values of f at time `t + Deltat` using the value `((delf)/(delt))` we found.
`((delf)/(delt))_(i,j)^(t)`
`overlinef_(i,j)^(t+Deltat)=f_(i,j)^(t)+ ((delf)/(delt))_(i,j)^(t)Deltat`
This value we found is called predicted value of f
Corrector step: In predictorstep we find the value of `(overline(delf)/(delt))`at time `t+Deltat` by rearward difference
we than find the average and use it in equation 1
`((delf)/(delt))_(avg) = 1/2[((delf)/(delt))_(i,j)^(t)+(overline(delf)/(delt))_(i,j)^(t+Deltat)]`
Boundary conditions:
Non-conservative form:
Nozzle shape and initial condition:
`A = 1 + 2.2(x-1.5)^2 0 le x le 3`
`T = 1 = 0.2314x`
`rho = 1 - 0.3146x`
`V = (0.1 + 1.09x)T^(1/2)`
Boundary conditions:
subsonic inflow(inlet):
`rho_1 = 1`
`T_1`=1
Supersonic outflow boundary (Ponit N):
`V_N = 2V_(N-1) - V_(N-2)`
`rho_N = 2rho_(N-1) - rho_(N-2)`
`V_N = 2V_(N-1) - V_(N-2)`
Conservative Form :
initial conditions:
Boudary conditions(Subsonic inlet):
`U_(1(i=1))=(rho^'A^')_(i=1)= (A^')_(i=1)`=fixed value
`U_(2(i=1))=2U_(2(i=2))-U_(2(i=3))`
`(U_3 = U_1((T^')/(gamma-1)+gamma/2V^'))_(i=1)`
supersonic outlet(N):
`(U_1)_N = 2(U_1)_(N-1) - (U_1)_(N-2)`
`(U_2)_N = 2(U_2)_(N-1) - (U_2)_(N-2)`
`(U_3)_N = 2(U_3)_(N-1) - (U_3)_(N-2)`
Time step for both methods: We calculate time step for every node at each iteration and choose the minimum among those.
`((Deltat)^t)_i = C (Deltax)/((a^t+V^t)_i`
At an adjacent grid point we have
`((Deltat)^t)_(i+1) = C (Deltax)/((a^t+V^t)_(i+1)`
Similarly we find the time step at each node and find the minimum among those:
`Deltat`=minimum`(Deltat_1,Deltat_2,Deltat_3.......Deltat_i....Deltat_N)^t`
Code:
We create functions for both conservative and non conservative form.
Non-conservative function:
function [mass_flow_rate,mass_flow_rate_t] = non_conservative(n)
%inputs
x = linspace(0,3,n);
dx = x(2)-x(1);
gamma = 1.4;
C = 0.5
nt = 1:1400;
%calculate initial profile
rho = 1-0.3146*x;
T = 1-0.2314*x;
v = (0.1+1.09*x).*T.^0.5;
A = 1+2.2*(x-1.5).^2;
M = v./(T.^0.5);
mass_flow_rate = rho.*A.*v;
p = rho.*T;
th = find(A==1)
rho_t = rho(th);
T_t = T(th);
v_t = v(th);
A_t = A(th);
M_t = M(th);
mass_flow_rate_t = mass_flow_rate(th);
p_t = p(th);
a= sqrt(T);
%time steps
nt = 1:1400;
for i=1:n
dt(i) = (C*dx)/(v(i)+a(i))
end
dt = min(dt);
%outer time loop
tic
for k = 1:length(nt)
%predictor method
rho_old = rho;
T_old = T;
v_old =v;
for j = 2:n-1
drhodx = (rho(j+1)-rho(j))/dx;
dvdx = (v(j+1)-v(j))/dx;
dTdx = (T(j+1)-T(j))/dx;
dlnAdx = (log(A(j+1))-log(A(j)))/dx;
%continuity equation
drhodt_p(j) = -rho(j)*dvdx-rho(j)*v(j)*dlnAdx-v(j)*drhodx;
%momentum equation
dvdt_p(j) = -v(j)*dvdx-(1/gamma)*(dTdx + (T(j)/rho(j))*drhodx);
%Energy equation
dTdt_p(j) = -v(j)*dTdx-(gamma-1)*T(j)*(dvdx+v(j)*dlnAdx);
rho(j) = rho_old(j)+drhodt_p(j)*dt;
v(j) = v_old(j)+drhodt_p(j)*dt;
T(j) = T_old(j)+dTdt_p(j)*dt;
%correction step
drhodx = (rho(j)-rho(j-1))/dx;
dvdx = (v(j)-v(j-1))/dx;
dTdx = (T(j)-T(j-1))/dx;
dlnAdx = (log(A(j))-log(A(j-1)))/dx;
%continuity equation
drhodt_c(j) = -rho(j)*dvdx-rho(j)*v(j)*dlnAdx-v(j)*drhodx;
%momentum equation
dvdt_c(j) = -v(j)*dvdx-(1/gamma)*(dTdx + (T(j)/rho(j))*drhodx);
%Energy equation
dTdt_c(j) = -v(j)*dTdx-(gamma-1)*T(j)*(dvdx+v(j)*dlnAdx);
%average
drhodt_avg(j)=0.5*(drhodt_c(j)+drhodt_p(j));
dvdt_avg(j)=0.5*(dvdt_c(j)+dvdt_p(j));
dTdt_avg(j)=0.5*(dTdt_c(j)+dTdt_p(j));
%updated parameter values
rho(j) = rho_old(j)+drhodt_avg(j)*dt;
v(j) = v_old(j)+dvdt_avg(j)*dt;
T(j) = T_old(j)+dTdt_avg(j)*dt;
end
%boundary conditions
v(1)=2*v(2)-v(3);
%outlet
v(31)=2*v(30)-v(29);
rho(31)=2*rho(30)-rho(29);
T(31)=2*T(30)-T(29);
p = rho.*T;
M = v./(T.^0.5);
mass_flow_rate = rho.*A.*v;
rho_t(k) = rho(th);
T_t(k) = T(th);
v_t(k) = v(th);
A_t(k) = A(th);
M_t(k) = M(th);
mass_flow_rate_t(k) = mass_flow_rate(th);
p_t(k) = p(th);
nt = 1:1400;
for i=1:n
dt(i) = (C*dx)/(v(i)+a(i))
end
dt = min(dt);
if n == 31
figure(1)
hold on
if k == 50
plot(x,mass_flow_rate,'b','LineWidth',1.75);
else if k == 100
plot(x,mass_flow_rate,'g','LineWidth',1.75);
else if k == 150
plot(x,mass_flow_rate,'r','LineWidth',1.75);
else if k == 200
plot(x,mass_flow_rate,'m','LineWidth',1.75);
else if k == 700
plot(x,mass_flow_rate,'-','LineWidth',1.75);
end
end
end
end
end
end
end
non_conservative_studytime = toc;
if n == 31
title('comparison of mass flow rate at different time steps')
xlabel('x/L (Non-dimensional distnace through th nozzle)')
ylabel('rhoVA/rho_{0}v_{0}A_{0}')
legend('50Deltat','100Deltat','150Deltat','200Deltat','700Deltat')
axis([0 3 0 2.5])
grid on
gridtext=sprintf('grid number %d',n)
saveas(gca,'gridtext_mass_time_non.jpeg')
hold off
end
if n == 31
figure(2)
yyaxis right
plot(x,M,'LineWidth',1.75);
title('variation of mach number and density along the nozzle')
xlabel('x/L (Non-dimensional distnace through th nozzle)')
ylabel('Mach number')
axis([0 3 0 3.5])
hold on
yyaxis left
plot(x,rho,'LineWidth',1.75);
xlabel('x/L (Non-dimensional distnace through the nozzle)')
ylabel('Nondimensional density rho/rho_{0}')
axis([0 3 0 1])
saveas(gca,'machdensity_non.jpeg')
hold off
figure(3)
subplot(4,1,1)
plot(nt,rho_t)
axis([0 1400 0 2])
ylabel('rho/rho_{0}')
subplot(4,1,2)
plot(nt,T_t)
axis([0 1400 0 2])
ylabel('T/T_{0}')
subplot(4,1,3)
plot(nt,p_t)
axis([0 1400 0 2])
ylabel('p/p_{0}')
subplot(4,1,4)
plot(nt,M_t)
axis([0 1400 0 2])
ylabel('M')
sgtitle('variation of primitive variables at throat at differnt time step')
saveas(gca,'nt_throat_non.jpeg')
figure(4)
subplot(4,1,1)
plot(x,rho)
axis([0 3 0 2])
ylabel('rho/rho_{0}')
subplot(4,1,2)
plot(x,T)
axis([0 3 0 2])
ylabel('T/T_{0}')
subplot(4,1,3)
plot(x,p)
axis([0 3 0 2])
ylabel('p/p_{0}')
subplot(4,1,4)
plot(x,M)
axis([0 3 0 4])
ylabel('M')
xlabel('x/L (Non-dimensional distnace through the nozzle)')
sgtitle('variation of primitive variables along the length')
saveas(gca,'x_primitive_non.jpeg')
end
end
Conservative method function:
function [m,m_t] = conservative(n)
L = 3;
x = linspace(0,L,n);
nt = 1:1400;
C = 0.5;
dx = L/(n-1);
g = 1.4;
rho = ones(1,n);
T = ones(1,n);
for i = 1:length(x)
if x(i) <= 0.5
rho(i) = 1;
T(i) = 1;
else if x(i) >= 0.5 && x(i)<=1.5
rho(i) = 1.0-0.366*(x(i)-0.5);
T(i) = 1.0-0.167*(x(i)-0.5);
else if x(i) >= 1.5 && x(i)<=3.5
rho(i) = 1.0-0.3879*(x(i)-1.5);
T(i) = 1.0-0.3507*(x(i)-1.5);
end
end
end
end
A = 1+(2.2*((x-1.5).^2));
%velocity
v = 0.59./(rho.*A);
%pressure
p = rho.*T
%mass flow rate
m = rho.*A.*v;
%Mach number
M = v./(sqrt(T));
th = find(A==1)
rho_t = rho(th);
T_t = T(th);
v_t = v(th);
A_t = A(th);
M_t = M(th);
m_t = m(th);
p_t = p(th);
a = sqrt(T);
for i=1:n
dt(i) = (C*dx)/(v(i)+a(i))
end
dt = min(dt);
u1 = rho.*A;
u2 = rho.*A.*v;
u3 = rho.*((T./(g-1)) + (g./2).*(v.^2)).*A;
tic
for k = 1:length(nt)
a = sqrt(T);
for i=1:n
dt(i) = (C*dx)/(v(i)+a(i));
end
dt = min(dt);
rho_old = rho;
T_old = T;
v_old =v;
f1 = u2;
f2 = ((u2.^2)./u1) + ((g-1)./g).*(u3 - (g/2).*((u2).^2./u1));
f3 = ((g.*u2.*u3)./u1) - ((g.*(g-1))./2).*((u2.^3)./(u1.^2));
%storing old values of u1,u2 and u3
u1_old = u1;
u2_old = u2;
u3_old = u3;
%predictor step
for i = 2:n-1
%source term
j2(i)=(1/g)*rho(i)*T(i)*((A(i+1)-A(i))/dx);
%continuity equation
du1dtp(i) = -(f1(i+1)-f1(i))./dx;
%momentum equation
du2dtp(i) = j2(i)-(f2(i+1)-f2(i))./dx;
%Energy equation
du3dtp(i)= -(f3(i+1)-f3(i))./dx;
%updated
u1(i) = (u1_old(i)+(du1dtp(i)*dt));
u2(i) = (u2_old(i)+(du2dtp(i)*dt));
u3(i) = (u3_old(i)+(du3dtp(i)*dt));
%predicted primitive variables
rho(i)=u1(i)./A(i);
v(i)=(u2(i)./u1(i));
T(i)=(g-1).*((u3(i)./u1(i)) - (g./2).*((u2(i)./u1(i)).^2));
p(i)=rho(i).*T(i);
end
%flux update
f1 = u2;
f2 = ((u2.^2)./u1) + ((g-1)./g).*(u3 - (g/2).*((u2).^2./u1));
f3 = ((g.*u2.*u3)./u1) - ((g.*(g-1))./2).*((u2.^3)./(u1.^2));
%corrector method
for i = 2:n-1
%source term
j2(i)=(1/g)*rho(i)*T(i)*((A(i)-A(i-1))/dx);
%continuity equation
du1dtc(i) = -(f1(i)-f1(i-1))./dx;
%momentum equation
du2dtc(i) = j2(i)-(f2(i)-f2(i-1))./dx;
%Energy equation
du3dtc(i)= -(f3(i)-f3(i-1))./dx;
end
%average values
du1dt_avg = 0.5*(du1dtc+du1dtp);
du2dt_avg = 0.5*(du2dtc+du2dtp);
du3dt_avg = 0.5*(du3dtc+du3dtp);
%solution update
for j = 2:n-1
u1(j)=(u1_old(j)+(du1dt_avg(j).*dt));
u2(j)=(u2_old(j)+(du2dt_avg(j).*dt));
u3(j)=(u3_old(j)+(du3dt_avg(j).*dt));
end
%finding values of primitive variables as they are updated for next
%iteration
rho = u1./A;
v =u2./u1;
T=(g-1).*((u3./u1) - (g./2).*((u2./u1).^2));
p=rho.*T;
m = rho.*A.*v;
M = v./(sqrt(T));
%applying boundary conditions
%inlet at node 1
u1(1) = rho(1)*A(1);
u2(1) = 2*u2(2) - u2(3);
u3(1) = u1(1)*((T(1)/(g-1)) + ((0.5*g)* (v(1)^2)));
%outlet at node n
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);
%flow variables at throat
rho_t(k) = rho(th);
T_t(k) = T(th);
v_t(k) = v(th);
M_t(k) = M(th);
m_t(k) = m(th);
p_t(k) = p(th);
%comparison of mass flow rate at different time steps
figure(6)
hold on
if k == 50
plot(x,m,'b','LineWidth',1.75);
else if k == 100
plot(x,m,'g','LineWidth',1.75);
else if k == 150
plot(x,m,'r','LineWidth',1.75);
else if k == 200
plot(x,m,'m','LineWidth',1.75);
else if k == 700
plot(x,m,'-','LineWidth',1.75);
end
end
end
end
end
end
conservative_study_time = toc;
title('comparison of mass flow rate at different time steps')
xlabel('x/L (Non-dimensional distnace through th nozzle)')
ylabel('rhoVA/rho_{0}v_{0}A_{0}')
legend('50Deltat','100Deltat','150Deltat','200Deltat','700Deltat')
axis([0 3 0 2.5])
grid on
saveas(gca,'gridtext_mass_time_conserv.jpeg')
hold off
figure(7)
yyaxis right
plot(x,M,'LineWidth',1.75);
xlabel('x/L (Non-dimensional distnace through th nozzle)')
ylabel('Mach number')
axis([0 3 0 3.5])
hold on
yyaxis left
plot(x,rho,'LineWidth',1.75);
xlabel('x/L (Non-dimensional distnace through th nozzle)')
ylabel('Nondimensional density rho/rho_{0}')
axis([0 3 0 1])
saveas(gca,'machdensity_conservative.jpeg')
hold off
figure(8)
title('variation of primitive variables at throat at differnt time step')
subplot(4,1,1)
plot(nt,rho_t)
axis([0 1400 0 2])
ylabel('rho/rho_{0}')
subplot(4,1,2)
plot(nt,T_t)
axis([0 1400 0 2])
ylabel('T/T_{0}')
subplot(4,1,3)
plot(nt,p_t)
axis([0 1400 0 2])
ylabel('p/p_{0}')
subplot(4,1,4)
plot(nt,M_t)
axis([0 1400 0 2])
ylabel('M')
sgtitle('variation of primitive variables at throat at differnt time step')
saveas(gca,'nt_throat_conserv.jpeg')
figure(9)
title('variation of primitive variables along the length of nozzle')
subplot(4,1,1)
plot(x,rho)
axis([0 3 0 2])
ylabel('rho/rho_{0}')
subplot(4,1,2)
plot(x,T)
axis([0 3 0 2])
ylabel('T/T_{0}')
subplot(4,1,3)
plot(x,p)
axis([0 3 0 2])
ylabel('p/p_{0}')
subplot(4,1,4)
plot(x,M)
axis([0 3 0 4])
ylabel('M')
xlabel('x/L (Non-dimensional distnace through the nozzle)')
sgtitle('variation of primitive variables along the length')
saveas(gca,'x_primitive_conservative.jpeg')
end
We now create main program and call two functions we have created.
clear all
close all
clc
%inputs
n = [31 61];
C = 0.5;
nt = 1:1400;
for i = 1:length(n)
x = linspace(0,3,n(i))
[mass_flow_rate_nonconservative,massflowrate_non_throat] = non_conservative(n(i))
[mass_flow_rate_conservative,massflowrate_conservative_throat] = conservative(n(i))
if n(i) == 31
figure(5)
subplot(2,1,1)
ylabel('rhoVA/rho_{0}v_{0}A_{0}')
xlabel('Time steps')
plot(nt,massflowrate_conservative_throat,'b',nt,massflowrate_non_throat,'r','LineWidth',1)
legend('conservative','Non conservative')
grid on
subplot(2,1,2)
plot(x,mass_flow_rate_conservative,'b',x,mass_flow_rate_nonconservative,'r','LineWidth',1)
ylabel('rhoVA/rho_{0}v_{0}A_{0}')
xlabel('length of the nozzle')
axis([0 3 0 1])
legend('conservative','Non conservative')
grid on
sgtitle('conservative vs nonconservative mass flow rate')
saveas(gcf,'nonvscon31.jpeg')
end
figure(10)
subplot(2,1,1)
title('Nonconservative mass flow rate vs time steps')
plot(nt,massflowrate_non_throat,'LineWidth',1)
ylabel('rhoVA/rho_{0}v_{0}A_{0}(Non-conservative)')
xlabel('Time steps')
hold on
grid on
legend('31 nodes','61 nodes')
subplot(2,1,2)
title('conservative mass flow rate vs time steps')
plot(nt,massflowrate_conservative_throat,'LineWidth',1)
ylabel('rhoVA/rho_{0}v_{0}A_{0}(conservative)')
xlabel('Time steps')
hold on
grid on
legend('31 nodes','61 nodes')
sgtitle('mass flow rate at throat')
saveas(gcf,'mvsnt.jpeg')
figure(11)
subplot(2,1,1)
plot(x,mass_flow_rate_nonconservative,'LineWidth',1)
ylabel('rhoVA/rho_{0}v_{0}A_{0}(non-conservative)')
xlabel('length of the nozzle')
axis([0 3 0 1])
hold on
legend('31 nodes','61 nodes')
subplot(2,1,2)
plot(x,mass_flow_rate_conservative,'LineWidth',1)
ylabel('rhoVA/rho_{0}v_{0}A_{0}(conservative)')
xlabel('Length of the nozzle')
axis([0 3 0 1])
hold on
legend('31 nodes','61 nodes')
sgtitle('mass flow rate along the length')
saveas(gcf,'mvsx.jpeg')
endResults:
Non-conservative form:
1.Variation of mass flow rate at throat for differnt time steps.
Here the nondimensional mass flow is plotted against nondimensional distance along the nozzle for different time steps. We can observe that 50 time steps the mass flow varies considerabaly and we can see that massflow rate is fluctuating inside nozzle this is due to transient variation of flow field variables. However after 200 time steps the mass flow distribution is beginning to settle down and after 700 time steps the mass fow is straight,horizontal across the graph. at this stage mass flow rate begins to converge
2. Variation of primitive variables along the length.
This shows the variation of primitive variables along the length of the nozzle as the flow moves from inlet to outlet.
3. Variation of primitive variables at the throat for differnt time steps.
It shows the variations of `rho`,p,T and M at the nozzle throat plotted versus the bumber of time steps. The abscissa starts at zero , which represents the initial conditions, and ends at time step 1400. We see that largest changes takes palce at early stages after which the steady state value is approached almost asymptotically.
4. Variation of mass flow rate and density along the nozzle.
conservative form:
1.Variation of mass flow rate at throat for differnt time steps.
2. Variation of primitive variables along the length.
3. Variation of primitive variables at the throat for differnt time steps.
4. Variation of mass flow rate and density along the nozzle.
Here dotted line indicates the adjusted density variation for increasing the node number from 61 to 31.
Conservative vs Non-conservative form:
First subplot : It shows the variation of mass flow rate at throat for different time steps for conservative and non conservative form, We can clearly see the conservative form is more stable and consistent from the start.
Second subplot: It shows the variation of mass flow rate along the length the nozzle for for conservative and non conservative form, both are stable but we can see little bumps in non conservative form.
Grid dependency test:
1. Mass flow rate at throat for differnt time steps for different number of nodes
2. Mass flow rate along the length for different number of nodes
conclusion :
We can clearly see that mass flow rate is almost same for differnt number of nodes therefore we can say that the methods are independent of number of grids and results always remain consistent.