Week 7 - Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
Here we consider the steady state isentrpic flow through convergent-divergent Nozzle. The flow at inlet is from reservoir. The flow expands isentropically to supersonic speeds at nozzle exit, where we have flow properties like Mach number,pressure,Temperature,velocity. The flow is locally subsonic at convergent section is nozzle ,sonic at throat, and supersonic at divergent section. The flow properties at Throat section is calculated and plot has been attached as well.
Even if the flow is 2D we make the assumption that flow properties vary along with x , this is tantamount to assuming uniform properties across any given section such a flow is Quasi-1D Flow.
To solve this we will set up a time-marching, finite-difference solution for the quasi 1D nozzle flow. The governing equation such as continuity,momentum,energy equations wrote in the term of PDE's suitable for time marching solution.
I solved this problem in both Non-conservation and Conservation form
The optimal criteria such as time steps, courant number, grids points, dt are taken.
1. In this projects i solved the 1D nozzle flow by implementing both non-conservation and conservation form of governing equations.
2. Performed Grid Dependence Test.
3. find Mass flow rate for both the terms.
4. Compared all the flow field properties using plot.
5. Distibution of non dimentional mass flow rate through the nozzle at different time steps.
Plot of Conservation form :
It shown that at initial time marching as the time derivatives are quite large the mass flow values appears to be very high through out the length of the nozzle. But after 200 time steps the value is approaching steady state and at 700 time steps we can say that the value is coverged. Mass flow rate does not change and remain constant between 700 to 1400 time steps.
Comparison the mass flow rate variations at different time steps
1. Mass flow rate distribution having less fluctuation in conservation form instead of non-conservation form.
2. The reason is that the mass flow rate is one of the main dependent variable in conservation form where as in non-conservation form dependent variable are primitives variable.
3. The conservation form generates more accurate results due to flux variable.
It seems that in figure 3 & 5 the forms return identical profile for mach number, denisty,temperature and pressure. The mass flow rate is more unstable in non conservation form.
In Figure 2 & 4 -
These graph are plotted for density,temperature,pressure and mach number.
Initially the values of the parameters are high and oscillate, the occurs due to various effects on flow such as expansion and compression which occurs in nozzle.
And in upcoming steps the values are approaching near to stability and contant condition.
It can be seen that there are more oscillations occurs in 0-200 times steps and graduallly decreases from 200-400 time steps.
From the plot we can also see that the number of cycles to run is 600 for convergance.
It can be said that the conservation form has less oscillation than non conservation form.
Conclusion -
1. Simulation time increase with increase in number of grid points but accuracy increases.
2. Non conseravtion form is faster than conservation form.
3. Mass flow distribution is better in conservation form than in non conseravation form as conservation form conserves more efficiency better.
4. Initially at lower time steps there are oscillations in the solution but as the time steps increases the solution moves to steady state and converges near to 800 time steps.
5. The simulation must run at least upto 800 time steps for it to converge.
Personally i imported lotof data from JD Anderson CFD basic book
% CONSERVATION FORM
function [m_c,M_c,v_c,T_c,P_c,rho_c,T_th_c,rho_th_c,P_th_c,M_th_c]=cons_solver(x,dx,n,C,g,nt)
%%
%Determining the initial profile of the density and T_c:
for i= 1:n
if(x(i)>=0 && x(i)<=0.5)
U1(i)=1;
U3(i)=1;
elseif(x(i)>=0.5 && x(i)<=1.5)
U1(i)=1-0.366*(x(i)-0.5);
U3(i)=1-0.167*(x(i)-0.5);
elseif(x(i)>=1.5 && x(i)<=3)
U1(i)=0.634-0.3879*(x(i)-1.5);
U3(i)=0.833-0.3507*(x(i)-1.5);
end
end
A=1+(2.2*(x-1.5).^2);%Area profile
th=find(A==1);%throat profile
U2=0.59./(U1.*A); %initial v_c profile
P_c=U1.*U3;%Initial profile of pressure ratio
a=sqrt(U3);%sonic v_c
M_c=U2./(U3.^0.5);%Mach number
m_c=U1.*U2.*A;%Mass flow rate
%Calculate the Initial solution vectors(U)
U1=U1.*A;
U2=U1.*A.*U2;
U3=U1.*A.*(U3./(g-1))+((g/2).*U2.^2);
%%
A=1+(2.2*(x-1.5).^2);
th=find(A==1);%throat profile
for i=1:n
if x(i)<=0.5
rho_c(i) = 1;
T_c(i) =1;
v_c(i) = 0.59/(rho_c(i)*A(i));
end
if x(i)>0.5 && x(i)<1.5
rho_c(i) = 1-(0.366*(x(i)-0.5));
T_c(i) = 1-(0.167*(x(i)-0.5));
v_c(i) = 0.59/(rho_c(i)*A(i));
end
if x(i)>=1.5
rho_c(i) = 0.634-(0.3879*(x(i)-1.5));
T_c(i) = 0.833-(0.3507*(x(i)-1.5));
v_c(i) = 0.59/(rho_c(i)*A(i));
end
end
a=sqrt(T_c);
U1 = rho_c.*A;
U2 = U1.*v_c;
U3 = U1.*((T_c./(g-1)) +((g/2)*v_c.^2));
v_c = v_c;
%outer time loop
for k=1:nt
dt=min(C.*dx./(a+v_c));%time step
U1_old=U1;
U2_old=U2;
U3_old=U3;
%Calculatingthe Flux vector(F)
F1=U2;
F2=((U2.^2)./U1)+((g-1)/g)*(U3-(g/2)*((U2.^2)./U1));
F3=((g*U2.*U3)./U1)-(g*(g-1)*(U2.^3)./(2*U1.^2));
%predictor step
for i=2:n-1
%Calculating the source terms:
J2=(1/g)*rho_c(i)*T_c(i)*(A(i+1)-A(i))/dx;
%continuity equation
dU1_dt_p(i)=-(F1(i+1)-F1(i))/dx;
%Momentum equation
dU2_dt_p(i)=(-(F2(i+1)-F2(i))/dx)+J2;
%Energy equation
dU3_dt_p(i)= -(F3(i+1)-F3(i))/dx;
%solution update
U1(i)=U1(i)+dU1_dt_p(i)*dt;
U2(i)=U2(i)+dU2_dt_p(i)*dt;
U3(i)=U3(i)+dU3_dt_p(i)*dt;
end
%Updating flux vectors
F1=U2;
F2=((U2.^2)./U1)+((g-1)/g)*(U3-(g/2)*((U2.^2)./U1));
F3=((g*U2.*U3)./U1)-(g*(g-1)*(U2.^3)./(2*U1.^2));
rho_c=U1./A;
v_c=U2./U1;
T_c=(g-1)*((U3./U1-(g/2)*(v_c).^2));
%corrector step
for i=2:n-1
%Updating the source terms:
J2=(1/g)*rho_c(i)*T_c(i)*(A(i)-A(i-1))/dx;
%continuity equation
dU1_dt_c(i)=-(F1(i)-F1(i-1))/dx;
%Momentum equation
dU2_dt_c(i)=-((F2(i)-F2(i-1))/dx)+J2;
%Energy equation
dU3_dt_c(i)= -(F3(i)-F3(i-1))/dx;
end
%Averaging the time derivatives
dU1_dt_avg=0.5*(dU1_dt_p+dU1_dt_c);
dU2_dt_avg=0.5*(dU2_dt_p+dU2_dt_c);
dU3_dt_avg=0.5*(dU3_dt_p+dU3_dt_c);
% final solution update
for i=2:n-1
U1(i)=U1_old(i)+dU1_dt_avg(i)*dt;
U2(i)=U2_old(i)+dU2_dt_avg(i)*dt;
U3(i)=U3_old(i)+dU3_dt_avg(i)*dt;
end
%
% %Final Primitive value update
% rho_c=U1./A;
% v_c=U2./U1;
% T_c=(g-1)*((U3./U1-(g/2)*(v_c).^2));
P_c=rho_c.*T_c;%pressure ratio
M_c=v_c./(sqrt(T_c));%Mach number
m_c=rho_c.*A.*v_c;%Mass flow rate
%Applying BC
%Inlet
U1(1)=rho_c(1)*A(1);
U2(1)=2*U2(2)-U2(3);
v_c(1)=U2(1)./U1(1);
U3(1)=U1(1)*((T_c(1)/(g-1))+((g/2)*(v_c(1)).^2));
%outlet
U2(n)=2*U2(n-1)-U2(n-2);
U1(n)=2*U1(n-1)-U1(n-2);
U3(n)=2*U3(n-1)-U3(n-2);
rho_c=U1./A;
v_c=U2./U1;
T_c=(g-1)*((U3./U1-(g/2)*(v_c).^2));
%initializing index for updation of primitive variable at throat wrt to
%increment in time step
rho_th_c(k)=rho_c(th);
T_th_c(k)=T_c(th);
P_th_c(k)=P_c(th);
M_th_c(k)=M_c(th);
v_th(k)=v_c(th);
m_th(k)=m_c(th);
%Comparision of non dimensional mass flow rates at different time:
mass_c=figure(7);
if k==50
plot(x,m_c,'linewidth',1.75);
hold on
elseif k==100
plot(x,m_c);
hold on
elseif k==200
plot(x,m_c);
hold on
elseif k==400
plot(x,m_c);
hold on
elseif k==700
plot(x,m_c);
hold on
elseif k==900
plot(x,m_c);
hold on
elseif k==1200
plot(x,m_c);
plot(x(end),m_c(end))
title('Comparision of Non-Dimensional mass flow rates at different time instances');
xlabel('x/L (Non-Dimensional distance through the nozzle');
ylabel('ND mass flow rate');
legend('50deltat','100deltat','200deltat','400deltat','700deltat','900deltat','1200deltat');
axis([0 3 0 2.5]);
grid on;
grid minor;
end
end
%
end # MAIN CODE
clear all
close all
clc
%Inputs
L=3; %Domain length
n=31;%no. grid points
th=(n+1)/2;% grid position of throat
x=linspace(0,L,n);
dx=x(2)-x(1);%mesh
%rho_stag=1; temp_stag=1;
g =1.4;%gamma
%M=0;
A=1+(2.2*(x-1.5).^2);%Area
%b=a*-1;
C=0.5;%CFL number
nt=1400;%no. of time steps
t_array=1:nt;
%solution for non conservative form:
tic;
[m,M,v,T,P,rho,T_th,rho_th,P_th,M_th]=NC_solver(x,dx,n,C,g,nt);
sim_time_NC=toc;
fprintf('The computation time for non conservative form =%f',sim_time_NC);
%Time wise variation of mach number density, pressure and temparature for Non-consrvative form
Time_variation=figure(2)
subplot(4,1,1)
plot(t_array,rho_th);
xlabel('Time')
ylabel('Density');
subplot(4,1,2)
plot(t_array,T_th);
xlabel('Time')
ylabel('Temperature');
subplot(4,1,3)
plot(t_array,P_th);
xlabel('Time')
ylabel('pressure');
subplot(4,1,4)
plot(t_array,M_th);
xlabel('Time')
ylabel('Mach number');
legend('density','temperature','pressure','Mach number');
%Qualitative aspect of 1D supersonic nozzle flow: Non conservative from:
Qualitative_C=figure(3)
subplot(4,1,1)
plot(x,rho);
xlabel('x')
ylabel('Density');
subplot(4,1,2)
plot(x,T);
xlabel('x')
ylabel('Temperature');
subplot(4,1,3)
plot(x,P);
xlabel('x')
ylabel('pressure');
subplot(4,1,4)
plot(x,M);
xlabel('x')
ylabel('Mach number');
%%
%solution for conservative form:
tic;
[m_c,M_c,v_c,T_c,P_c,rho_c,T_th_c,rho_th_c,P_th_c,M_th_c]=cons_solver(x,dx,n,C,g,nt);
sim_time_C=toc;
fprintf('The computation time for non conservative form =%f',sim_time_C);
%Time wise variation of mach number density, pressure and temparature for consrvative form
Time_variation=figure(4)
subplot(4,1,1)
plot(t_array,rho_th_c);
xlabel('Time')
ylabel('Density');
subplot(4,1,2)
plot(t_array,T_th_c);
xlabel('Time')
ylabel('Temperature');
subplot(4,1,3)
plot(t_array,P_th_c);
xlabel('Time')
ylabel('pressure');
subplot(4,1,4)
plot(t_array,M_th_c);
xlabel('Time')
ylabel('Mach number');
legend('density','temperature','pressure','Mach number');
%Qualitative aspect of 1D supersonic nozzle flow: conservative from:
Qualitative_C=figure(5)
subplot(4,1,1)
plot(x,rho_c);
xlabel('x')
ylabel('Density');
subplot(4,1,2)
plot(x,T_c);
xlabel('x')
ylabel('Temperature');
subplot(4,1,3)
plot(x,P_c);
xlabel('x')
ylabel('pressure');
subplot(4,1,4)
plot(x,M_c);
xlabel('x')
ylabel('Mach number');
# NON-CONSERVATION FORM
function [m,M,v,T,P,rho,T_th,rho_th,P_th,M_th]=NC_solver(x,dx,n,C,g,nt)
%Initial conditions(at t=0)
rho_i=1-0.3146*x; %Density
T_i=1-0.2314*x;%Temperature
v_i=(0.1+1.09*x).*(T_i.^0.5);%velocity
A=1+(2.2*(x-1.5).^2);%Area
%nd_x=x/L
%Fielf variables
rho=rho_i;T=T_i;v=v_i;
a=sqrt(T);%sonic velocity
P=rho.*T;%pressure ratio
M=v./(T.^0.5);%Mach number
m=rho.*v.*A;%Mass flow rate
%stencil metrics for throat section
th=find(A==1);
%calculating time step
% for i=1:n
% dt_calc(i)=(C*dx)/(a(i)+v(i));
% end
% dt=min(dt_calc);
%outer time loop
for k=1:nt
dt=min(C.*dx./(a+v));
%initial conditions as old conditions
rho_old=rho; T_old=T;v_old=v;
%predictor step
for i=2:n-1
drho_dx=(rho(i+1)-rho(i))/dx;
dv_dx=(v(i+1)-v(i))/dx;
dT_dx=(T(i+1)-T(i))/dx;
dlogA_dx= (log(A(i+1))-log(A(i)))/dx;
%continuity equation
drho_dt_p(i)=-(rho(i)*dv_dx) - (rho(i)*v(i)*dlogA_dx) - (v(i)*drho_dx);
%Momentum equation
dv_dt_p(i)=-(v(i)*dv_dx) - ((1/g)*(dT_dx + (T(i)*drho_dx/rho(i))));
%Energy equation
dT_dt_p(i)= -(v(i)*dT_dx) - ((g-1)*T(i)*(dv_dx + (v(i)*dlogA_dx)));
%solution update
rho(i)=rho(i)+drho_dt_p(i)*dt;
v(i)=v(i)+dv_dt_p(i)*dt;
T(i)=T(i)+dT_dt_p(i)*dt;
end
%corrector step
for i=2:n-1
drho_dx=(rho(i)-rho(i-1))/dx;
dv_dx=(v(i)-v(i-1))/dx;
dT_dx=(T(i)-T(i-1))/dx;
dlogA_dx= (log(A(i))-log(A(i-1)))/dx;
%continuity equation
drho_dt_c(i)=-(rho(i)*dv_dx) - (rho(i)*v(i)*dlogA_dx) - (v(i)*drho_dx);
%Momentum equation
dv_dt_c(i)=-(v(i)*dv_dx) - ((1/g)*(dT_dx + (T(i)*drho_dx/rho(i))));
%Energy equation
dT_dt_c(i)= -(v(i)*dT_dx) - ((g-1)*T(i)*(dv_dx + (v(i)*dlogA_dx)));
end
%Averaging the time derivatives
drho_dt_avg=0.5*(drho_dt_p+drho_dt_c);
dv_dt_avg=0.5*(dv_dt_p+dv_dt_c);
dT_dt_avg=0.5*(dT_dt_p+dT_dt_c);
% final solution update
for i=2:n-1
rho(i)=rho_old(i)+drho_dt_avg(i)*dt;
v(i)=v_old(i)+dv_dt_avg(i)*dt;
T(i)=T_old(i)+dT_dt_avg(i)*dt;
end
%Applying BC
%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);
P=rho.*T;%pressure ratio
a=sqrt(T);
M=v./a;%Mach number
m=rho.*v.*A;%Mass flow rate
%initializing index for updation of primitive variable at throat wrt to
%increment in time step
rho_th(k)=rho(th);
T_th(k)=T(th);
P_th(k)=P(th);
M_th(k)=M(th);
v_th(k)=v(th);
m_th(k)=m(th);
%Comparision of non dimensional mass flow rates at different time:
mass_nc=figure(1);
if k==50
plot(x,m,'linewidth',1.75);
hold on
elseif k==100
plot(x,m);
hold on
elseif k==200
plot(x,m);
hold on
elseif k==400
plot(x,m);
hold on
elseif k==700
plot(x,m);
hold on
elseif k==900
plot(x,m);
hold on
elseif k==1200
plot(x,m);
plot(x(end),m(end))
title('Comparision of Non-Dimensional mass flow rates at different time instances');
xlabel('x/L (Non-Dimensional distance through the nozzle');
ylabel('ND mass flow rate');
legend('50deltat','100deltat','200deltat','400deltat','700deltat','900deltat','1200deltat');
axis([0 3 0 2.5]);
grid on;
grid minor;
end
end