Week 7 - Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
The quasi one dimensional nozzle flow Mackormack 1D flow is considered only along the x direction. The nozzle is a convergent divergent nozzle. The mach number at the throat
Equations for the non conservative form of 1D Mackomarck method problem:
continuity equation:
`dρ/dt=−ρ∗dv/dx−ρ∗v∗d(lnA)/dx−v∗dρ/dx`
momentum equation:
`dv/dt=−v∗dv/dx−1/γ∗(dt/dx+t∗dρ/dx/ρ)dv/dt`
energy equation:
`dT/dt=−v∗dT/dx−(γ−1)∗T∗(dv/dx+v∗d(lnA)/dx)dT/dt`
For the non conservative form of the equation:
U1 = ρ∗A
U2 = ρ∗A∗v
U3 =ρ∗(T/(γ−1)+γ/2∗v2)∗A
F1 =ρ∗A∗v
`F2 = ρ∗A∗v2+1∗P∗A/(γ)\rho*A*v^2`
`F3 = ρ∗(T/(γ−1)+γ∗v2/2)∗V∗A+P∗A∗V`
`J2 = 1/(γ)∗P∗dA/dx`
`dU1/dt=−dF1/dx`
`dU2/dt=−dF2/dx + J2`
ρ,P,T,V are primitive terms. The flux terms F1,F2,F3 could directly be expressed in terms of these primitive variables
but the solution would blow up in such a case in the subsonic region after an estimated of 300 iterations and
as a result the flux terms are expressed in the form of source terms U1, U2,U3. If the flux terms are expressed
in the form of primitive terms just like written above, there is a lack of purity in the solution and the solution
would blow up .The below equations represent F1,F2,F3 completely in terms of U1,U2,U3.
F1 = U2
`F2 = (U2^2)/(U1)+((γ−1)/γ)∗(U3−0.5∗γ∗(U2^2)/(U1))`
n=61;
x= linspace(0,3,n);
dx= x(2)-x(1);
gamma= 1.4;
c=0.5;
nt=1400;
a_con= 1+ 2.2*(x-1.5).^2;
a= 1+ 2.2*(x-1.5).^2;
tic
% function calling
[rho, v, t,mfr, p, mach, rho_throat, v_throat, t_throat, p_throat, mfr_throat, mach_throat]= nonconservative(n,x,gamma,nt,c);
timenon= toc;
tic
%function calling
[throat,rho_con, v_con, t_con,p_con,p_con_throat,v_con_throat,rho_con_throat,mfr_con_throat,t_con_throat,mach_con_throat, mfr_con, mach_con]= conservative(x,c,nt,gamma,n);
timecon= toc;
figure(3)
subplot(2,3,1)
plot(x,rho,'r', x, rho_con, 'b*');
xlabel('non dimensional nozzle length');
ylabel('non dimensional density');
grid minor
legend('non conservative','conservative');
title('density plots');
subplot(2,3,2)
plot(x,v,'r', x, v_con, 'b*');
xlabel('non dimensional nozzle length');
ylabel('non dimensional velocity');
grid minor
legend('non conservative','conservative');
title('velocity plots');
subplot(2,3,3)
plot(x,t,'r', x, t_con, 'b*');
xlabel('non dimensional nozzle length');
ylabel('non dimensional temperature');
grid minor
legend('non conservative','conservative');
title('temperature plots');
subplot(2,3,4)
plot(x,p,'r', x, p_con, 'b*');
xlabel('non dimensional nozzle length');
ylabel('non dimensional pressure');
grid minor
legend('non conservative','conservative');
title('pressure plots');
subplot(2,3,5)
plot(x,mfr,'r', x, mfr_con, 'b*');
xlabel('non dimensional nozzle length');
ylabel('non dimensional mass flow rates');
grid minor
legend('non conservative','conservative');
title('mass flow rate plots');
subplot(2,3,6)
plot(x,mach,'r', x, mach_con, 'b*');
xlabel('non dimensional nozzle length');
ylabel('non dimensional mach number');
grid minor
legend('non conservative','conservative');
title('mach number plots');
figure(4)
subplot(2,3,1)
plot(1:nt,rho_throat,'r',1:nt, rho_con_throat,'b');
xlabel('time step ');
ylabel('throat density');
grid minor
legend('non conservative','conservative');
title('throat density plots');
axis([0 1400 0 1]);
%figure(5)
subplot(2,3,2)
plot(1:nt,v_throat,'r',1:nt, v_con_throat,'b');
xlabel('time step ');
ylabel('throat velocity');
grid minor
legend('non conservative','conservative');
title('throat velocity plots');
axis([0 1400 0.7 1.4]);
%figure(6)
subplot(2,3,3)
plot(1:nt,t_throat,'r',1:nt, t_con_throat,'b');
xlabel('time step ');
ylabel('throat temperature');
grid minor
legend('non conservative','conservative');
title('throat temperature plots');
axis([0 1400 0.66 0.96]);
%figure(7)
subplot(2,3,4)
plot(1:nt,p_throat,'r',1:nt, p_con_throat,'b');
xlabel('time step ');
ylabel('throat pressure');
grid minor
legend('non conservative','conservative');
title('throat pressure plots');
axis([0 1400 0.3 1]);
%figure(8)
subplot(2,3,5)
plot(1:nt,mfr_throat,'r',1:nt, mfr_con_throat,'b');
xlabel('time step ');
ylabel('throat mass flow rate');
grid minor
legend('non conservative','conservative');
title('throat mass flow rate plots');
axis([0 1400 0.4 1.1]);
%figure(9)
subplot(2,3,6)
plot(1:nt,mach_throat,'r',1:nt, mach_con_throat,'b');
xlabel('time step ');
ylabel('throat mach no');
grid minor
legend('non conservative','conservative');
title('throat mach no plots');
axis([0 1400 0.8 1.8]);
function [rho, v, t,mfr, p, mach, rho_throat, v_throat, t_throat, p_throat, mfr_throat, mach_throat]= nonconservative(n,x,gamma,nt,c)
n=61;
x= linspace(0,3,n);
dx= x(2)-x(1);
% calculate initial profiles
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;
gamma= 1.4;
c=0.5;
for xyz= 1:n
if x(xyz)==1.5
throat= xyz;
end
end
% time steps
dt = min((c.*dx)./((t.^0.5) + v));
nt= 1400;
%outer time loop
for k= 1:nt
rho_old= rho;
t_old= t;
v_old= v;
%predictor method
for j= 2:n-1
dvdx= (v(j+1)-v(j))/dx;
drhodx= (rho(j+1)-rho(j))/dx;
dtdx= (t(j+1)-t(j))/dx;
dlogadx= (log(a(j+1))-log(a(j)))/dx;
%predicting the continuity equation of the non conservative form
drhodt_p(j) = -rho(j)*dvdx -rho(j)*v(j)*dlogadx - v(j)*drhodx;
%predicting the momentum equation of the non conservative form
dvdt_p(j)= -v(j)*dvdx -(1/gamma)*(dtdx + (t(j)/rho(j))*drhodx);
%predicting the energy equation of the non conservative form
dtdt_p(j)= -v(j)*dtdx -(gamma-1)*t(j)*(dvdx + v(j)*dlogadx);
% prediction is done in forward differencing
% solution update is done inside this loop as we can do element by
% element updating at each of the inner nodes.
rho(j)= rho(j)+ drhodt_p(j)*dt;
v(j)= v(j)+ dvdt_p(j)*dt;
t(j)= t(j)+ dtdt_p(j)*dt;
end
for j= 2:n-1
% corrector method
dvdx= (v(j)-v(j-1))/dx;
drhodx= (rho(j)-rho(j-1))/dx;
dtdx= (t(j)-t(j-1))/dx;
dlogadx= (log(a(j))-log(a(j-1)))/dx;
%correcting the continuity equation of the non conservative form
drhodt_c(j) = -rho(j)*dvdx -rho(j)*v(j)*dlogadx - v(j)*drhodx;
%correcting the momentum equation of the non conservative form
dvdt_c(j)= -v(j)*dvdx -(1/gamma)*(dtdx + (t(j)/rho(j))*drhodx);
%correcting the energy equation of the non conservative form
dtdt_c(j)= -v(j)*dtdx -(gamma-1)*t(j)*(dvdx + v(j)*dlogadx);
% correction is done in backward differencing
end
% average values
drhodt= 0.5*(drhodt_p + drhodt_c);
dvdt= 0.5*(dvdt_p + dvdt_c);
dtdt= 0.5*(dtdt_p + dtdt_c);
% final answers
for i= 2:n-1
rho(i)= rho_old(i) + drhodt(i)*dt;
v(i)= v_old(i) + dvdt(i)*dt;
t(i)= t_old(i) + dtdt(i)*dt;
end
% applying boundarybconditions
v(1)= 2*v(2)-v(3); % inlet
rho(1)=1;
t(1)=1;
v(n)= 2*v(n-1)- v(n-2); %outlet
rho(n)= 2*rho(n-1)- rho(n-2);
t(n)= 2*t(n-1)- t(n-2);
% calculating the pressure mach number and mass flow rate
mfr= rho.*v.*a; % mass flowrate
p= rho.*t; % pressure
mach= v./t.^0.5;
% calculating the throat values after every time step
rho_throat(k)= rho(throat);
v_throat(k)= v(throat);
t_throat(k)= t(throat);
p_throat(k)= p(throat);
mfr_throat(k)= mfr(throat);
mach_throat(k)= mach(throat);
% plotting the mass flowrates at different time steps
b= figure(1)
hold on;
if k==1
plot(x, mfr, 'b');
end
if k==50
plot(x, mfr,'r');
end
if k ==100
plot(x, mfr, 'g');
end
if k== 250
plot(x, mfr,'y');
end
if k==500
plot(x,mfr,'p');
end
if k==1000
plot(x,mfr,'c');
end
legend('time step 1','time step 50','time step 100','time step 250','time step 500','time step 1000');
xlabel('non dimensional length x/l');
ylabel('non dimensional mass flow rate');
title('non conservative form mass flow rates');
end
endfunction [throat,rho_con, v_con, t_con,p_con,p_con_throat,v_con_throat,rho_con_throat,mfr_con_throat,t_con_throat,mach_con_throat, mfr_con, mach_con]= conservative(x,c,nt,gamma,n)
clear all
close all
clc
gamma= 1.4;
n=61;
c=0.5;
x= linspace(0,3,n);
dx= x(2)-x(1);
%rho_con= ones(1,n);
%t_con= ones(1,n);
a_con= 1+ 2.2*(x-1.5).^2;
% defining the rho and t
for xy= 1:n
if (x(xy)>=0 && x(xy)<0.5)
rho_con(xy)=1;
t_con(xy)=1;
end
if (x(xy)>=0.5 && x(xy)< 1.5)
rho_con(xy) = 1- 0.366*(x(xy) - 0.5);
t_con(xy)= 1- 0.167*(x(xy) - 0.5);
end
if (x(xy)>=1.5 && x(xy)< 3.5)
rho_con(xy) = 0.634- 0.3879*(x(xy) - 1.5);
t_con(xy)= 0.833- 0.3507*(x(xy) - 1.5);
end
end
v_con= 0.59./(rho_con.*a_con) ;
u1 = rho_con.*a_con;
u2= rho_con.*a_con.*v_con;
u3= rho_con.*a_con.*((t_con)/(gamma-1) + ((gamma/2).*v_con.^2));
% finding the time step
nt= 1400; % number of time steps
dt = min((c.*dx)./((t_con.^0.5) + (v_con)));
for th= 1:n
if x(th)== 1.5
throat= th;
end
end
for pq= 1:nt
u1_old= u1;
u2_old= u2; % preserving the old values of u
u3_old = u3;
% expressing f1,f2,f3 and j2 in terms of solution vector or pure
% form
f1_p= u2;
f2_p= (u2.^2./u1) + (gamma-1)/(gamma)*(u3- (0.5*gamma.*u2.^2)./u1);
f3_p= (gamma*(u2.*u3))./u1 - (0.5*gamma*(gamma-1).*u2.^3)./(u1.^2);
% for the inner nodes now we start the predictor step
for j= 2:n-1
df1pdx= (f1_p(j+1)- f1_p(j))/dx;
df2pdx= (f2_p(j+1)- f2_p(j))/dx;
df3pdx= (f3_p(j+1)- f3_p(j))/dx;
du1dt_p(j) = -df1pdx; % continnuit equation
du2dt_p(j) = -df2pdx + (1/(gamma))*rho_con(j)*t_con(j)*((a_con(j+1)-a_con(j))/dx); % p= rho_t., momentum equation
du3dt_p(j)= -df3pdx; % energy equation
% updating the predicted values
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 of the rho, v,t,p values after predictor step
rho_con= u1./a_con;
v_con= u2./u1;
t_con= (gamma-1)*((u3./u1)- (gamma/2).*v_con.^2);
p_con= rho_con.*t_con;
% updating the values of f1, f2, f3 that is to be used for
% corrector step
f1_c= u2;
f2_c= (u2.^2./u1) + (gamma-1)/(gamma).*(u3- (0.5*gamma.*u2.^2)./u1);
f3_c= (gamma.*(u2.*u3))./u1 - (0.5*gamma*(gamma-1).*u2.^3)./(u1.^2);
% corrector step backward differencing
for ab= 2:n-1
df1cdx= (f1_c(ab)- f1_c(ab-1))/dx;
df2cdx= (f2_c(ab)- f2_c(ab-1))/dx;
df3cdx= (f3_c(ab)- f3_c(ab-1))/dx;
du1dt_c(ab) = -df1cdx; % continnuity equation
du2dt_c(ab) = -df2cdx + (1/(gamma))*rho_con(ab)*t_con(ab)*((a_con(ab)-a_con(ab-1))/dx); % p= rho_t., momentum equation
du3dt_c(ab)= -df3cdx; % energy equation
end
% average values of dudts
du1dt_avg= 0.5.*(du1dt_p + du1dt_c);
du2dt_avg= 0.5.*(du2dt_p + du2dt_c);
du3dt_avg= 0.5.*(du3dt_p + du3dt_c);
% updating the final values after the corrector step
for k= 2:n-1
u1(k)= u1_old(k)+ du1dt_avg(k)*dt;
u2(k)= u2_old(k)+ du2dt_avg(k)*dt;
u3(k)= u3_old(k)+ du3dt_avg(k)*dt;
end
% inlet boundary conditions at the subsonic inlet
u1(1)= rho_con(1)*a_con(1);
u2(1)= 2*u2(2)- u2(3);
v_con(1)= u2(1)/u1(1);
u3(1)= u1(1)*(t_con(1)/(gamma-1) + 0.5*gamma*v_con(1)^2);
%outlet boundary conditions
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);
% update of the rho, v,t,p values after corrector step
rho_con= u1./a_con;
v_con= u2./u1;
t_con= (gamma-1).*((u3./u1)- (gamma/2).*v_con.^2);
% mass flow rate, pressure and mach number
p_con= rho_con.*t_con;
mfr_con = rho_con.*a_con.*v_con;
mach_con= v_con./t_con.^0.5;
% mass flow rate, p and mach number at the throat for each timestep
p_con_throat(pq)= p_con(throat);
v_con_throat(pq)= v_con(throat);
t_con_throat(pq)= t_con(throat);
rho_con_throat(pq)= rho_con(throat);
mfr_con_throat(pq)= mfr_con(throat);
mach_con_throat(pq)= mach_con(throat);
%plotting the mass flowrates at different time steps
a= figure(2)
hold on;
if pq==1
plot(x, mfr_con, 'b');
end
if pq==50
plot(x, mfr_con,'r');
end
if pq ==100
plot(x, mfr_con, 'g');
end
if pq== 250
plot(x, mfr_con,'y');
end
if pq==500
plot(x,mfr_con,'p');
end
if pq==1000
plot(x,mfr_con,'c');
end
legend('time step 1','time step 50','time step 100','time step 250','time step 500','time step 1000');
xlabel('non dimensional length x/l');
ylabel('non dimensional mass flow rate');
title('conservative form mass flow rates');
axis([0 3 0.3 1.1]);
end
end Grid Independence test for n=61:
Results of grid independence studies for n= 31 and n=61 for both the conservation and non conservation forms
| Density ratio | temperature ratio | pressure ratio | velocity ratio | mass flow rate ratio | |
| n=31, conservative form |
0.6524
|
0.8410
|
0.5487
|
0.8969
|
0.5851
|
| n=31, non-conservative form |
0.6387
|
0.8365
|
0.5342
|
0.9140
|
0.5838
|
| n=61, conservative form | 0.9690 | 0.9875 | 0.9568 | 0.2688 | 0.5829 |
| n=61, non-conservative form | 0.9691 | 0.9875 | 0.9570 | 0.2690 | 0.5832 |
Conclusions
The above table shows the values of the various variables at the point i=16 for both the conservative and non conservative forms. It indicates that for number of nodes 31 and 61, the results obtained are the same in case of both conservative and non conservative form. The grid independence is valid and the solution doesnt change with the mesh.
Both the conservative and non conservative forms give similar results and both the forms of the equations are one and the same. It is as per the problem and the problem specification, the appropriate form has to be used.
At the throat of the nozzle, for both the grids, it can be observed that the values of each of the parameters(density, mass flow rate,pressure,temperature,velocity and mach number), almost achieve the state of convergence and their values do not change much after a certain point of time. The curve gets straight and almost parallel to the x axis. This means that the solution attained steady state.
Even though the solution gets converged in almost 1000 iterations in the above graphs, it is good to calculate the solution for atleast 1400 to 1600 iterations as proper amount of convergence with conclusive results could further be obtained