Week 7 - Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
The objective is to solve steady, isentropic flow through a convergent-divergent nozzle as shown:
The flow at the inlet to the nozzle comes from a reservoir where the pressure and temperature are denoted by p0 and T0, respectively. It is assumed that at a given section,where the area is A, the flow properties remain constant. This type of flow is defined as quasi-one-dimensional flow.
The governing equations can be derived in two forms, conservative and non-conservative. The equations obtained from a finite volume fixed in space in either integral or differential form, are called the conservation form of governing equations. The equations obtained from a finite control volume moving with the fluid, are called as non-conservation form of governing equation.
The governing equations for non-conservative form for the quasi-1D flow as described above are as follows:
Continuity equation:
`(∂(ρA))/(∂t)+(∂(ρAV))/(∂x)=0`
Momentum Equation:
`ρA (∂V)/(∂t)+ρAV (∂V)/(∂x) =-A (∂ρ)/(∂x)`
Energy Equation:
`ρC_v (∂T)/(∂t)+ρVC_v (∂T)/(∂x)=-ρRT*((∂V)/(∂x)+V (∂(lnA))/(∂x))`
For nozzle flows, the flow-field variables are frequently expressed in terms of non-dimensional variables, where the flow referenced to their reservoir values. The non-dimensional variables vary between 0 and 1, proving to be beneficial when presenting results.
So the non-dimensional are given as:
`ρ^'=ρ/ρ_0`
`T^'=T/T_0`
Dimensional length is given as:
`x^'=x/L`
Denoting speed of sound in reservoir as `a_0=√(γRT_o )`
Dimensionless velocity is given as:
`v^'=v/a_0`
Also dimensionless time is given as:`t^'=t/(L/a_0 )`
And dimensionless area: `A^'=A/A^*`
Where A* is the sonic throat area.
So the governing equations in the form of the above dimensionless variable become,
Continuity:`(∂ρ^')/(∂t^' )=-ρ^' (∂V^')/(∂x^' )-ρ^' V^' (∂(lnA^'))/(∂x^' )-V'(∂ρ^')/(∂x')`
Momentum:`(∂V^')/(∂t^' )=-V^' (∂V^')/(∂x^' )-(1/γ)*((∂T^')/(∂x^' )+T^'/ρ^'(∂ρ^')/(∂x^' ))`
Energy:`(∂T^')/(∂t^' )=-V^' (∂T^')/(∂x^' )-(γ-1)T'((∂V^')/(∂x^' )+V'(∂(lnA^'))/(∂x^' ))`
The above equations are then solved through discretizing them in finite difference forms.
The MacCormack’s technique is used to solve the equations using time marching approach. The technique consists of predictor step using forward differencing for spatial derivatives and then predicting values of the variables. These predicted values are then used in the corrector step which uses backward differencing for spatial derivatives. The advantage of using this technique is that it achieves second order accuracy without the hassle of calculating the terms up to second order derivative in Taylor series, which usually involve more algebra.
So for the predictor step:
`((∂ρ)/(∂t))_i^t=-ρ_i^t (V_(i+1)^t-V_i^t)/(Δx)-ρ_i^t V_i^t (lnA_(i+1)-lnA_i)/(Δx)-V_i^t (ρ_(i+1)^t-ρ_i^t)/(Δx)`
`((∂V)/(∂t))_i^t=-V_i^t (V_(i+1)^t-V_i^t)/(Δx)-(1/γ)*((T_(i+1)^t-T_i^t)/(Δx)+(T_i^t)/(ρ_i^t ) (ρ_(i+1)^t-ρ_i^t)/(Δx))`
`((∂T)/(∂t))_i^t=-V_i^t (T_(i+1)^t-T_i^t)/(Δx)-(γ-1) T_i^t ((V_(i+1)^t-V_i^t)/(Δx)+V_i^t (lnA_(i+1)-lnA_i)/(Δx))`
We can then obtain predicted values as:
`ρ ̅_i^(t+Δt)=ρ_i^t+((∂ρ)/(∂t))_i^t Δt`
`V ̅_i^(t+Δt)=V_i^t+((∂V)/(∂t))_i^t Δt`
`T ̅_i^(t+Δt)=T_i^t+((∂T)/(∂t))_i^t Δt`
Then for corrector step
`((∂rho ̅)/(∂t))_i^(t+Δt)=-ρ ̅_i^(t+Δt) (V ̅_i^(t+Δt)-V ̅_(i-1)^(t+Δt))/(Δx)-ρ ̅_i^(t+Δt) V ̅_i^(t+Δt) (lnA_i-lnA_(i-1))/(Δx)-V ̅_i^(t+Δt) (ρ ̅_i^(t+Δt)-ρ ̅_(i-1)^(t+Δt))/(Δx)`
`((∂V ̅)/(∂t))_i^(t+Δt)=-V ̅_i^(t+Δt) (V ̅_i^(t+Δt)-V ̅_(i-1)^(t+Δt))/(Δx)-(1/γ)*((T ̅_i^(t+Δt)-T ̅_(i-1)^(t+Δt))/(Δx)+(T ̅_i^(t+Δt))/(ρ ̅_i^(t+Δt) ) (ρ ̅_i^(t+Δt)-ρ ̅_(i-1)^(t+Δt))/(Δx))`
`((∂T ̅)/(∂t))_i^(t+Δt)=-V ̅_i^(t+Δt) (T ̅_i^(t+Δt)-T ̅_(i-1)^(t+Δt))/(Δx)-(γ-1) T ̅_i^(t+Δt) ((V ̅_i^(t+Δt)-V ̅_(i-1)^(t+Δt))/(Δx)+V ̅_i^(t+Δt) (lnA_i-lnA_(i-1))/(Δx))`
Now the average time derivative is given as:
`((∂ρ)/(∂t))_(av)=0.5[((∂ρ)/(∂t))_i^t+((∂rho ̅)/(∂t))_i^(t+Δt) ]`
`((∂V)/(∂t))_(av)=0.5[((∂V)/(∂t))_i^t+((∂V ̅)/(∂t))_i^(t+Δt) ]`
`((∂T)/(∂t))_(av)=0.5[((∂T)/(∂t))_i^t+((∂T ̅)/(∂t))_i^(t+Δt) ]`
Finally, the corrected values,
`ρ_i^(t+Δt)=ρ_i^t+((∂ρ)/(∂t))_(av) Δt`
`V_i^(t+Δt)=V_i^t+((∂V)/(∂t))_(av) Δt`
`T_i^(t+Δt)=T_i^t+((∂T)/(∂t))_(av) Δt`
The time step is calculated according to the Courant number criterion:
`Δt=C (Δx)/(a+V)`
Where C is the courant number `<=`1
Even though the system which is being solved is non linear, the courant number criterion does provide good guidance. In the given case the vaue of courant number is taken to be 0.5.
As ‘a’ which is speed of sound and V keep on changing throughout the grid, the timestep come out to be different for each grid point. The least value of these is considered to be the final value of time step.
Boundary Conditions:
Point 1 as shown in the figure previously, is essentially in the reservoir. It represents a subsonic inflow boundary condition with flow coming from the reservoir and entering the nozzle. In exit, point N is an outflow boundary, with flow leaving the nozzle at exit.
So for subsonic inflow (point 1), the boundary conditions are as follows,
`V_1=2V_2-V_3`
`ρ_1=1`
`T_1=1`
So for supersonic outflow (point N), the boundary conditions are as follows,
`V_N=2V_(N-1)-V_(N-2)`
`ρ_N=2ρ_(N-1)-ρ_(N-2)`
`T_N=2T_(N-1)-T_(N-2)`
Initial Conditions:
The area of nozzle varies with respect to x as:
`A=1+2.2(x-1.5)^2` `0≤x≤3`
`ρ=1-0.3146x`
`T=1-0.2314x`
`V=(0.1+1.09x) T^(1/2)`
The results obtained are as follows:
In the above plot, the values for non dimensional mass flow are plotted for various timesteps across the length of nozzle. It is seen that the mass flow is almost a horizontal straight line after 700 timesteps which agrees with the fact about nozzles that the mass flow rate across nozzle is constant.
In the above graph, the steady state distribution of non dimensional quantities are plotted across the length of nozzle. It is seen that the steady state was reached after 477 timesteps for a tolerance criteria of 1e-4. Also the mach number is seen to be greater than 1 after the throat section which suggests that the flow is supersonic at the outlet of the nozzle.
Conservative Form:
The governing equations for conservative form for the quasi-1D flow in non-dimensional form are as follows:
Continuity:`(∂(ρ'A'))/(∂t')+(∂(ρ^' A'V'))/(∂x')=0`
Momentum:`(∂(ρ'A'V^'))/(∂t^' )+(∂[ρ'A'V'^2+(1/γ)ρ'A'])/(∂x^' )=(1/γ) ρ'(∂A')/(∂x')`
Energy:`(∂[ρ^' (e^'/(γ-1)+γ/2 V'^2 )A'])/(∂t^' )+(∂[ρ^' (e^'/(γ-1)+γ/2 V'^2 ) V^' A^'+ρ'A'V'])/(∂x^' )=0`
Now, we can define the elements of the solutions vector U, the flux vector F, the source term J as follows:
`U_1=ρ'A'`
`U_2=ρ'A'V'`
`U_3=ρ^' (e^'/(γ-1)+γ/2 V'^2 )A'`
`F_1=ρ^' A^' V^'=U_2`
`F_2=ρ'A'V'^2+(1/γ) ρ^' A^'=(U_2^2)/U_1 +((γ-1)/γ)[U_3-(γ/2) (U_2/U_1 )^2 ]`
`F_3=ρ^' (e^'/(γ-1)+γ/2 V'^2 ) V^' A^'+ρ^' A^' V^'=γ (U_2 U_3)/U_1 -(γ((γ-1)/2) U_2^3)/(U_1^2 )`
`J_2=(1/γ) ρ'(∂A')/(∂x')`
With these elements the governing equations become:
`(∂U_1)/(∂t')=-(∂F_1)/(∂x')`
`(∂U_2)/(∂t')=-(∂F_2)/(∂x^' )+J_2`
`(∂U_3)/(∂t')=-(∂F_3)/(∂x')`
The above equations are then discretized by Maccormack’s method by following same procedure as for non-conservative part.
The primitive variables are then calculated as:
`ρ^'=U_1/A^'`
`V^'=U_2/U_1`
`T^'=(γ-1)(U_3/U_1 -γ/2 V'^2 )`
`p^'=ρ'T'`
Boundary conditions:
Subsonic inflow:
`U_(1(i=1)) =A_(i=1)^'=fixed value`
`U_(2(i=1)) =2U_(2(i=2)) -U_(2(i=3) )`
`U_(3(i=1)) =U_1 (T^'/(γ-1)+γ/2 V'^2 )`
Supersonic outflow:
`U_(1N)=2U_(1(N-1))-U_(1(N-2))`
`U_(2N)=2U_(2(N-1))-U_(2(N-2))`
`U_(3N)=2U_(3(N-1))-U_(3(N-2))`
Initial Conditions:
`ρ^'=1` for `0≤x≤0.5`
`T^'=1`
`ρ^'=1-0.366(x^'-0.5)` for `0.5≤x≤1.5`
`T^'=1-0.167(x^'-0.5)`
`ρ^'=0.634-0.3879(x^'-1.5)` for `1.5≤x≤3.5`
`T^'=0.833-0.3507(x^'-1.5)`
`V^'=(U_2)/(ρ^' A^' )=0.59/(ρ'A')`
For time step, the same procedure is used as for non-conservative part.
The results obtained are as follows:
In the above plot, the values for non dimensional mass flow are plotted for various timesteps across the length of nozzle. It is seen that the mass flow is almost a horizontal straight line after 700 timesteps which agrees with the fact about nozzles that the mass flow rate across nozzle is constant for non conservative part
The above plots show the steady state distributions for the various primitive variables. The steady state is reached in 795 timesteps.
From the number of timesteps required, it can be concluded that the non conservative formulation was faster to converge to steady state.
Grid dependency Check
It was seen that the values changed substantially when making the grid finer and finer for the non conservative formulation. Thus the solution is not grid independent in case of non conservative part and the grid must be made finer and finer until grid convergence is achieved.
Unlike the non conservative form, the conservative form changed values very little on making the grid finer. For the practical purposes, it can be assumed that the conservative form was grid independent. However, as the grid is made finer the residuals drop very slowly for the conservative part.
Comparison of non dimesional mass flow:
In the above graph, the mass flow rates obtained from non conservative and the conservative form are plotted together. It can be seen that on a magnified scale, there are substantial fluctuations at the inlet and outlet of the non conservative part. However the conservative part shows almost a straight line which is representative of the actual flow inside a nozzle which indeed has constant mass flow rate.
In conclusion, for the given flow problem, we cannot establish a certain superiority of one form over another. The following observations can be made:
1. The conservation form yields a better mass flow distribution. It does a better job of conserving mass.
2. The non conservation form leads to smaller residuals and performs faster than the conservative part.
3. The calculation effort for the conservative part is certainly greater due to the need to decode the primitive and the flux variable; such decoding is not necessary when using non conservative part.
main program to call the functions:
% main program to call the functions
%Gaurav Gokhale
clear all
close all
clc
%inputs
n=31; %grid points
L=3; %Length
nt=1400; %timesteps
[rhonc,Tnc,vnc,Mnc,mfrnc]=nonconservative(n,L,nt);
[rhoc,Tc,vc,Mc,mfrc]=conservative(n,L,nt);
x=linspace(0,3,31);
figure(5)
hold on
plot(x,mfrc,'r')
plot(x,mfrnc,'b')
legend('conservation form','non conservation form')
title('Comparison of steady state Mass flow rate ')
ylabel('Non dimensional mass flow')
Funtion for non conservative form:
function [rho,T,v,M,mfr]=nonconservative(n,L,nt)
gamma=1.4;
C=0.5;
x=linspace(0,L,n);
dx=x(2)-x(1);
%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;
dt=zeros(1,length(a));
for m=1:length(a)
dt(m)=C*(dx/(sqrt(T(m))+v(m)));
end
dt=min(dt);
for k=1:nt
k;
rho_old=rho;
v_old=v;
T_old=T;
%Predictor step
for i=2:n-1
dvdx=(v(i+1)-v(i))/dx;
dTdx=(T(i+1)-T(i))/dx;
drhodx=(rho(i+1)-rho(i))/dx;
dadx=(log(a(i+1))-log(a(i)))/dx;
%continuity eqn
drhodtp(i)=-rho(i)*dvdx-rho(i)*v(i)*dadx-v(i)*drhodx;
%momentum eqn
dvdtp(i)=-v(i)*dvdx-(1/gamma)*(dTdx+(T(i)/rho(i))*drhodx);
%energy eqn
dTdtp(i)=-v(i)*dTdx-(gamma-1)*T(i)*(dvdx+v(i)*dadx);
%predicted values
rho(i)=rho_old(i)+drhodtp(i)*dt;
v(i)=v_old(i)+dvdtp(i)*dt;
T(i)=T_old(i)+dTdtp(i)*dt;
end
%Corrector step
for j=2:n-1
dvdx=(v(j)-v(j-1))/dx;
dTdx=(T(j)-T(j-1))/dx;
drhodx=(rho(j)-rho(j-1))/dx;
dadx=(log(a(j))-log(a(j-1)))/dx;
%continuity eqn
drhodtc(j)=-rho(j)*dvdx-rho(j)*v(j)*dadx-v(j)*drhodx;
%momentum eqn
dvdtc(j)=-v(j)*dvdx-(1/gamma)*(dTdx+(T(j)/rho(j))*drhodx);
%energy eqn
dTdtc(j)=-v(j)*dTdx-(gamma-1)*T(j)*(dvdx+v(j)*dadx);
end
%avg derivative
drhodt=0.5*(drhodtp+drhodtc);
dvdt=0.5*(dvdtp+dvdtc);
dTdt=0.5*(dTdtp+dTdtc);
%Corrected Values
for l=2:n-1
rho(l)=rho_old(l)+drhodt(l)*dt;
v(l)=v_old(l)+dvdt(l)*dt;
T(l)=T_old(l)+dTdt(l)*dt;
end
%Boundary Conditions
%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;
M=v./T.^(1/2);
mfr=rho.*a.*v;
figure(3)
suptitle(sprintf('No.of timesteps: %d',k))
% hold on
subplot(2,2,1)
plot(x,p,'y')
ylabel('Pressure')
axis([0 3 0 2.2])
title('Pressure Vs X')
subplot(2,2,2)
plot(x,T,'b')
ylabel('Temperature')
axis([0 3 0 2.2])
title('Temperature Vs X')
subplot(2,2,3)
plot(x,M,'r')
ylabel('Mach Number')
axis([0 3 0 4])
title('Mach Number Vs X')
subplot(2,2,4)
plot(x,rho,'g')
ylabel('Density')
axis([0 3 0 1.5])
title('Density Vs X')
% figure(4)
% hold on
% if k==50
% plot(x,mfr,'r')
% elseif k==100
% plot(x,mfr,'b')
% elseif k==150
% plot(x,mfr,'g')
% elseif k==200
% plot(x,mfr,'c')
% elseif k==700
% plot(x,mfr,'y')
% legend('50dt','100dt','150dt','200dt','700dt')
% title('Non dimensional mass flow at different timesteps non-conservative form')
% xlabel('Non dimensional distance through nozzle (x)')
% ylabel('Non Dimensional mass flow')
% end
% pause(0.03)
errorv=max(abs(v_old-v));
errort=max(abs(T_old-T));
errorrho=max(abs(rho_old-rho));
error=max([errorv,errort,errorrho]);
if error<1e-4
break
end
end
end
Function for conservative part:
function [rho,T,v,M,mfr]=conservative(n,L,nt)
%inputs
x=linspace(0,3,n);
dx=x(2)-x(1);
gamma=1.4;
C=0.5;
%initial Profiles
a=1+2.2*(x-1.5).^2;
for i=1:n
if x(i)<0.5
rho(i)=1;
T(i)=1;
elseif 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);
elseif 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
v(i)=0.59/(rho(i)*a(i));
end
U1=rho.*a;
U2=rho.*a.*v;
U3=rho.*((T./(gamma-1))+(gamma/2).*v.^2).*a;
% nt=1400;
dt=zeros(1,length(a));
for m=1:length(a)
dt(m)=C*(dx/(sqrt(T(m))+v(m)));
end
dt=min(dt);
% dt=0.0267;
for k=1:nt
k
U1_old=U1;
U2_old=U2;
U3_old=U3;
F1=U2;
F2=U2.^2./U1+((gamma-1)/gamma).*(U3-(gamma/2)*U2.^2./U1);
F3=gamma.*U2.*U3./U1-0.5*gamma*(gamma-1).*U2.^3./U1.^2;
%Predictor Step
for i=2:n-1
J2=(1/gamma)*rho(i)*T(i)*((a(i+1)-a(i))/dx);
dU1dtp(i)=-(F1(i+1)-F1(i))/dx;
dU2dtp(i)=-(F2(i+1)-F2(i))/dx+J2;
dU3dtp(i)=-(F3(i+1)-F3(i))/dx;
%predicted values
U1(i)=U1_old(i)+dU1dtp(i)*dt;
U2(i)=U2_old(i)+dU2dtp(i)*dt;
U3(i)=U3_old(i)+dU3dtp(i)*dt;
end
F1=U2;
F2=U2.^2./U1+((gamma-1)/gamma).*(U3-(gamma/2)*U2.^2./U1);
F3=gamma.*U2.*U3./U1-0.5*gamma*(gamma-1).*U2.^3./U1.^2;
rho=U1./a;
T=(gamma-1).*(U3./U1-(gamma/2).*(U2./U1).^2);
%Corrector Step
for i=2:n-1
J2=(1/gamma)*rho(i)*T(i)*((a(i)-a(i-1))/dx);
dU1dtc(i)=-(F1(i)-F1(i-1))/dx;
dU2dtc(i)=-(F2(i)-F2(i-1))/dx+J2;
dU3dtc(i)=-(F3(i)-F3(i-1))/dx;
end
%avg Derivatives
dU1dt=0.5*(dU1dtp+dU1dtc);
dU2dt=0.5*(dU2dtp+dU2dtc);
dU3dt=0.5*(dU3dtp+dU3dtc);
%Corrected Values
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
%Boundary conditions
%inlet
U1(1)=a(1);
U2(1)=2*U2(2)-U2(3);
v(1)=U2(1)/U1(1);
U3(1)=U1(1)*(T(1)/(gamma-1)+(gamma/2)*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);
%primitives
rho=U1./a;
v=U2./U1;
T=(gamma-1).*(U3./U1-(gamma/2).*v.^2);
M=v./T.^(1/2);
mfr=U2;
p=rho.*T;
figure(1)
suptitle(sprintf('No.of timesteps: %d',k))
% hold on
subplot(2,2,1)
plot(x,p,'y')
ylabel('Pressure')
axis([0 3 0 2.2])
title('Pressure Vs X')
subplot(2,2,2)
plot(x,T,'b')
ylabel('Temperature')
axis([0 3 0 2.2])
title('Temperature Vs X')
subplot(2,2,3)
plot(x,M,'r')
ylabel('Mach Number')
axis([0 3 0 4])
title('Mach Number Vs X')
subplot(2,2,4)
plot(x,rho,'g')
ylabel('Density')
axis([0 3 0 1.5])
title('Density Vs X')
% figure(2)
% hold on
% if k==50
% plot(x,U2,'r')
% elseif k==100
% plot(x,U2,'b')
% elseif k==150
% plot(x,U2,'g')
% elseif k==200
% plot(x,U2,'c')
% elseif k==700
% plot(x,U2,'y')
% legend('50dt','100dt','150dt','200dt','700dt')
% title('Non dimensional mass flow at different timesteps conservative form')
% xlabel('Non dimensional distance through nozzle (x)')
% ylabel('Non Dimensional mass flow')
% end
pause(0.03)
errorU1=max(abs(U1_old-U1));
errorU2=max(abs(U2_old-U2));
errorU3=max(abs(U3_old-U3));
error=max([errorU1,errorU2,errorU3])
if error<1e-4
break
end
end
end