Week 7 - Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method

The steady isoentropic flow through the convergent -divergent nozzle as given :

flow at the inlet to the nozzle comes from the reservoir is large (therotically `A ->oo`0 and hence the velocity is very small `(V->0)` .Thus ,`p_o` and `T_o` are the stagnantation values or total pressure and total temperature,respectively.The flow expands isoentropically to supersonic speeds at the nozzle exit,where the exit pressure,temperature where exit parameters are reperesented with subscript e.

The flow is locally subsonic in the convergent section of the nozzle,sonic at the throat (minimum area) ,and supersonic at the divergent section . The sonic flow (M=1)at the throat means the local velocity at this section is equal to the local speed of the sound.We can repersent the throat values with our own representation either with * or any other.

Some of the useful results obtained.

`(A/A^**)^2=1/M^2*(2/(gamma+1)*(1+(gamma-1)/2*M^2))^((gamma+1)/(gamma-1))`

where `gamma` is the ratio of specific heat.For air at standard conditons it is equal to 1.4.

Some other equations include

`(p/p_o)=(1+(gamma-1)/2*M^2)^(-gamma/(gamma-1)`

`(rho/rho_o)=(1+(gamma-1)/2*M^2)^(-1/(gamma-1)`

`(T/T_o)=(1+(gamma-1)/2*M^2)^(-1)`

Mach number:

Mach number is a dimensionless quantity in fluid dynamics representing the ratio of flow velocity past a boundry to the local speed of sound

`M=u/c ,`
where:

M is the local Mach number,
u is the local flow velocity with respect to the boundaries (either internal, such as an object immersed in the flow, or external, like a channel), and
c is the speed of sound in the medium, which in air varies with the square root of the thermodynamics temeprature.

Assuming air to be an ideal gas, the formula to compute Mach number in a subsonic compressible flow is derived from Bernoullis equation for M < 1:

`M=sqrt(2/(gamma+1)*[(q_c/p+1)^((gamma-1)/(gamma))-1]`

and the speed of sound varies with teperature as follows:

`c=sqrt(gamma*R_**T)`

where`q_c` is imact pressure

`p` is the static pressure.

`gamma` is the ratio of specific heat of a gas at a constant pressure to heat at a constant volume (1.4 of air).

`R_**` is the specific gas constant.

Set Up

MaCcormack mathod:

In computational fluid dynamics, the MacCormack method is a widely used discretization scheme for the numerical solution of hyperbolic partial differential equations.

Maccormacks method provides second order accuracy without the invlovement of second order derivative thus reducing the complexity of the calculation.It is a explicit method.

According to Maccormacks method the value of variable at grid point i,j at time t+dt is given as follows:

`F^(actual)= F^( Initial) + (del(f))/(delt)_(avg)*delt`

The average time derivative can be found as follows ,which is the half of sum of predicotor value and correcotr value .

`((delv)/(delt))_(avg)=((delv)/(delt))_(p)+((delv)/(delt))_(c)`

So basically before that we need to set up equations for solving the predictor and corrector method.

Governing equations:

The super sonic flow through nozzle is governed by Navier stokes equations,which could be derived for both the cases of conservative and Non conservative methods.

The continuity equation can be given as follows:

`(delrho^')/(delt^')=-(rho^')*(delV^')/(delx^')-rho^'*V^'*(del(lnA^'))/(delx')-V^'*((delrho^')/(delx^'))`

The momentum equation:

`(delV^')/(delt)^'=-V^'*(delV^')/(delx^')-1/gamma*((delT^')/(delx^')+T^'/(rho^')*(delrho^')/(delx^'))`

The energy equation is as follows:

`(delT^')/(delt^')=-V^'*(delT^')/(delx^')-(gamma-1)*T^'*[(delV^')/(delx^')+V^'*(del(lnA^'))/(delx^')]`

These equation can be derived from

The `T_o` and `rho_o` arethe reservoir temperature and density.

We define the non dimensional temperature and density respectively `

define the initial condition :

• The shape of nozzle ,A=A( x) is specified and held fixed ,independent of time,For this case we take the parabolic area distribution
• `A=1+2.2*(x-1.5)^2` `0<=x<=3`
• For density ,Temperature and velocity is given by
• `rho = 1-0.3146*x`
• `T=1-0.2314x`
• `V=(0.1+1.09x)*sqrt(T)`
• Pressure `P=rho*T`
• Mach number= `u/sqrtT`
• The above mentioned are the initial conditions at t=0

predictor step:

This is a kind of descritisation step.Where we descretize with forward difference equation.

We first set up the finite difference expressions as follows:

We set up the spatial derivatives as forward differences .Also to reduce the complexity of the notation we will drop the prime to denote the dimensionless variable.

`rid generation we divide the length of the nozzle in number of points along the axis .The characteristics of the flow changes from subsonic to super sonic from right to left .The distance between the two grid points is dx.

Smaller the mesh size the error also decreases`

In discretisation part we apply thr Macromach technique to solve the the conservative and non conservative forms of governing equations.

we get the predicted equation as follows: These are denoted by barred quantities

Equations `rho_i^t,V_i^t` and `T_i^t` are the known values at time t

The corrcetor method:

the backward differencing scheme and rest remains same as that of the predictor method.

Continuity equation:

Momentum equation:

Energy equation:

Calculation of time step:

dt =C*(dx)/(u+a)`

where C is the courant number itts value is 0.5

u is the local velocity

a is the local speed of sound

Boundry condtions at subsonic section

linear extrapolation

`V_1=2*V_2-V_3`

Density and Temperature at point 1 is held fixed.

Supersonic outflow boundry conditions:

`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)`

The non dimensional conservation form of the coninuity ,momentum and energy equations for quasi-one -dimensional flow ,respectively.The equations for quasi one dimensional flow can be expressed in generic form.Let us define the elements of the solution vector U ,flux vector F and the source term J as follows:

These are the equations to solve numerically by using the Macomacks technique.

To obtain the primitve variables we must decode the elements U1,U2,U3 as follows.From the defination of U1,U2 and U3 given above

Applying the boundary conditions:

Inlet boundry conditions &

Outlet boundry condition:The flow properties at the downstream ,supersonic outflow boundry are obtained by linear extrapolation from the two adjacent internal points.

If N denotes the grid point at the outflow boundary

values of F1,F2 and F3 at grid point i=N are obtained from the values of U1,U2 and U3 at point i=N using the equations.

Initial conditions:

Commands:

close all
clear all
clc

%input parameters
n=31;
x=linspace(0,3,n);
dx=x(2)-x(1);
gamma=1.4;
CFL=0.5;
Area=1+2.2*(x-1.5).^2;
throat=((n+1)/2);
nt=1400;

%main body
tic;
[V,rho,T,m,M,P,V_t,rho_t,T_t,m_t,M_t,P_t]=non_con(throat,x,gamma,dx,Area,n,nt,CFL);
non_con=toc;

tic;
[Vc,rhoc,Tc,cons_mass_flow_throat,Mc,pc,Vc_throat,rho_throat,T_throat,cons_mass_flow,Mc_throat,pc_throat]=con_ser(throat,x,gamma,dx,Area,n,nt,CFL);
con_ser=toc;

%plotting section
figure(3)
plot(x,m,'linewidth',0.95);
hold on;
plot(x,cons_mass_flow,'*');
%exact analytical solution rho*A*v=0.579
plot(x,linspace(0.579,0.579,n),'--','linewidth',1.5);
xlabel('Non Dimensional distance of nozzle (x)');
ylabel('mass flow rate');
legend('non conservative form','conservative form','exact analytical solution');
title('Comparision of mass flow rate(conventional form vs non conventional');
grid on ;

figure(4)
plot(x,M,'linewidth',1.5);
hold on
plot(x,Mc,'*','linewidth',0.65,'color','g');
grid on
xlabel('distance of nozzle(x)');
ylabel('Mach number m');
legend('Non conservative form','conservative form');
title('Mach number profile');
grid on ;

figure(5)
plot(x,P,'linewidth',1.25);
hold on
plot(x,pc,'g','linewidth',1.25);
grid on
xlabel('distance of nozzle(x)');
ylabel('pressure');
legend('non conservative form','conservative form');
title('Pressure Profile');
grid on;

figure(6)
plot(x,rho,'linewidth',1.5)
hold on

plot(x,rhoc,'*','linewidth',1.25);
grid on
xlabel('distance of nozzle(x)');
ylabel('density');
legend('non conservative form','conservative form');
title('density profile');
grid on;

figure(7)
plot(x,T,'linewidth',1.5);
hold on
plot(x,Tc,'*','linewidth',0.65,'color','g');
grid on
xlabel('distance of nozzle(x) ');
ylabel('Temperature');
legend('non conservative form','conservative form');
title('Temperature change');

figure(8)
plot(x,V,'linewidth',1.0);
hold on ;
plot(x,Vc,'*','linewidth',0.65,'color','g');
grid on
xlabel('distance of nozzle(x)');
ylabel('velocity');
legend('non conservative form','conservative form');
title('Velocity profile');
grid on;

figure(9)
plot(1:nt,rho_throat,'linewidth',1.5);
hold on
plot(1:nt,rho_t,'c','linewidth',1.65);
grid on
ylabel('density');
legend('Conservative form','non conservative form');
title('density variation at various time step');
grid on;

figure(10)
plot(1:nt,T_throat,'linewidth',0.75);
hold on
plot(1:nt,T_t,'linewidth',0.95)
grid on
ylabel('temeperature');
legend('Conservative form','non conservative form');
title('Temperature variation at various time step at throat');
grid on;

figure(11)
plot(1:nt,pc_throat,'linewidth',0.85)
hold on
plot(1:nt,P_t,'linewidth',0.95)
grid on
ylabel('Presssure');
legend('conservative form','non conservative form');
title('pressure variatiuon at various time step at throat');
grid on;

figure(12)
plot(1:nt,Mc_throat,'linewidth',0.85);
hold on
plot(1:nt,M_t,'linewidth',0.90)
grid on;

xlabel('Time step');
ylabel('Mach number');
legend('Conservative form','non conservative form');
title('Mach number variation at various time step at throat');

function command :

function for non-conservative form

function [V,rho,T,m,M,P,V_t,rho_t,T_t,m_t,M_t,P_t]=non_con(throat,x,gamma,dx,A,n,nt,C)

%calculations

rho =1-0.3146*x;
T=1-0.2314*x;
V=(0.1+1.09*x).*T.^0.5;

dt=min(C.*(dx./(sqrt(T)+V)));

for it=1:nt
rho_old=rho;
V_old=V;
T_old=T;

for j=2:n-1

drhodx_p=(rho(j+1)-rho(j))/dx;
dvdx_p=(V(j+1)-V(j))/dx;
dtdx_p=(T(j+1)-T(j))/dx;
dlogA_dx_p=(log(A(j+1))-log(A(j)))/dx;

%continuity equation
drhodt_p(j)=-rho(j)*dvdx_p-rho(j)*V(j)*dlogA_dx_p-V(j)*drhodx_p;
%momentum equation
dvdt_p(j)=-V(j)*dvdx_p-(1/gamma)*(dtdx_p+(T(j)/rho(j))*drhodx_p);
%energy equation
dtdt_p(j)=-V(j)*dtdx_p-(gamma-1)*T(j)*(dvdx_p+V(j)*dlogA_dx_p);
%solution update
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

%corrector step.
for j=2:n-1
drhodx_c=(rho(j)-rho(j-1))/dx;
dvdx_c=(V(j)-V(j-1))/dx;
dtdx_c=(T(j)-T(j-1))/dx;
dlogA_dx_c=(log(A(j))-log(A(j-1)))/dx;

%continuity equation
drhodt_c(j)=-(rho(j)*dvdx_c)-rho(j)*V(j)*dlogA_dx_c-V(j)*drhodx_c;
%momentum equation
dvdt_c(j)=-(V(j)*dvdx_c)-(1/gamma)*(dtdx_c+(T(j)/rho(j))*drhodx_c);
%energy equation
dtdt_c(j)=-(V(j)*dtdx_c)-(gamma-1)*T(j)*(dvdx_c+V(j)*dlogA_dx_c);
end

%average values
dvdt=0.5*(dvdt_c+dvdt_p);
drhodt=0.5*(drhodt_p+drhodt_c);
dtdt=0.5*(dtdt_c + dtdt_p);

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
V(1)=2*V(2)-V(3);
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);

m=rho.*A.*V;
P=rho.*T;
M=V./sqrt(T);


rho_t(it)=rho(throat);
V_t(it)=V(throat);
T_t(it)=T(throat);
P_t(it)=P(throat);
M_t(it)=M(throat);
m_t(it)=m(throat);

rho_exact=0.634*ones(1,1400);
t_exact=0.833*ones(1,1400);
p_exact = 0.528*ones(1,1400);
mach_exact=1*ones(1,1400);

figure(1)
if it==50
plot(x,m,'b','linewidth',1.5);
hold on
elseif it==100
plot(x,m,'r','linewidth',1.5);
hold on
elseif it==150
plot(x,m,'c','linewidth',1.5);
hold on
elseif it==300
plot(x,m,'g','linewidth',1.5);
hold on
elseif it==600
plot(x,m,'y','linewidth',1.5);
hold on
elseif it==900
plot(x,m,'k','linewidth',1.5);
hold on
elseif it==1400
plot(x,m,'m','linewidth',1.5);
grid on
end
end

xlabel('Non dimensional distance through nozzle (x)','fontsize',10);
ylabel('Non dimensional mass flow rate through nozzle x','fontsize',10);
title('Variation of mass flow rate (non conseravtive form) along x at different steps');
legend({'50 time steps','100 time steps','150 time steps','300 time steps','600 time steps','900 time steps','1400 time steps'});


end

function command for

Conservative function:

%conservative form
function [Vc,rhoc,Tc,cons_mass_flow_throat,Mc,pc,Vc_throat,rho_throat,T_throat,cons_mass_flow,Mc_throat,pc_throat]=con_ser(throat,x,gamma,dx,A,n,nt,C)

R=8.314;
for i=1:n
if (x(i)<=0.5)
rho_in (i)=1;
T_in(i)=1;
elseif (x(i)>=0.5) && (x(i) <= 1.5)
rho_in(i)=1-(0.366*(x(i)-0.5));
T_in(i)=1-0.167*(x(i)-0.5);
else
rho_in(i)=0.634-0.3879*(x(i)-1.5);
T_in(i)=0.833-0.3507*(x(i)-1.5);
end
end

u_in=0.59./(rho_in.*A);
rhoc = rho_in;
Tc=T_in;
Vc=u_in;


dt=0.0206;
t_sim=(nt-1)*dt;
t=0:dt:t_sim;

U1=rhoc.*A;
U2=rhoc.*A.*Vc;
U3=rhoc.*A.*((Tc/(gamma-1))+(0.5*gamma.*Vc.*Vc));

%Macormack method
for j=1:length(t)
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)-((gamma*(gamma-1))/2)*(U2.^3)./(U1.^2);

%predictor method
for k=2:n-1
dA_dx(j)=(A(k+1)-A(k))/dx;
dlogA_dx = (log(A(k+1))-log(A(k)))/dx;
J2(k)=dlogA_dx*((gamma-1)/gamma)*(U3(k)-(0.5*gamma*U2(k)^2)/(U1(k)));


%continuity step

dF1_dx_p=(F1(k+1)-F1(k))/dx;
dU1_dt_p(k)=-dF1_dx_p;
%momentum equation
dF2_dx_p=(F2(k+1)-F2(k))/dx;
dU2_dt_p(k)=-dF2_dx_p+J2(k);

dF3_dx_p=(F3(k+1)-F3(k))/dx;
dU3_dt_p(k)=-dF3_dx_p;

U1(k)=(dU1_dt_p(k).*dt)+U1_old(k);
U2(k)=(dU2_dt_p(k).*dt)+U2_old(k);
U3(k) = (dU3_dt_p(k).*dt)+U3_old(k);
end
%applying boundry condition
%input
U1(1)=A(1);
U2(1)=2*U2(2)-U2(3);
U3(1)=U1(1)*((Tc(1)/(gamma-1)+((gamma/2)*Vc(1)^2)));
%output
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);

rhoc=U1./A;
Tc=(gamma-1).*((U3./U1)-((gamma/2).*(U2./U1).^2));

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);

%corrector method
for k=2:n-1
dA_dx(j)=(A(k)-A(k-1))/dx;
dlogA_dx=(log(A(k))-log(A(k-1)))/dx;
J2(k)=dlogA_dx*((gamma-1)/gamma)*(U3(k)-(0.5*gamma*U2(k)^2)/(U1(k)));

%continuity step
dF1_dx_c=(F1(k)-F1(k-1))/dx;
dU1_dt_c(k)=-dF1_dx_c;

%momentum equation
dF2_dx_c=(F2(k)-F2(k-1))/dx;
dU2_dt_c(k)=-dF2_dx_c +J2(k);

%energy equation
dF3_dx_c=(F3(k)-F3(k-1))/dx;
dU3_dt_c(k)=-dF3_dx_c;
end
for pc=2:n-1
dU1_dt_avg(pc)=0.5*(dU1_dt_c(pc)+dU1_dt_p(pc));
dU2_dt_avg(pc)=0.5*(dU2_dt_c(pc)+dU2_dt_p(pc));
dU3_dt_avg(pc)=0.5*(dU3_dt_c(pc)+dU3_dt_p(pc));

U1(pc)=(dU1_dt_avg(pc)*dt)+U1_old(pc);
U2(pc)=(dU2_dt_avg(pc)*dt)+U2_old(pc);
U3(pc)=(dU3_dt_avg(pc)*dt)+U3_old(pc);
end

%boundry
%input
U1(1)=A(1);
U2(1)=2*U2(2)-U2(3);
U3(1)=U1(1)*((Tc(1)/(gamma-1)+((gamma/2)*Vc(1)^2)));
%input
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);


rhoc=U1./A;
Vc=U2./U1;
Tc=(gamma-1)*((U3./U1)-((gamma/2).*(Vc.^2)));

pc=rhoc.*Tc;
Mc=Vc./(Tc.^0.5);
cons_mass_flow=rhoc.*Vc.*A;
%variables at throat
rho_throat(j)=rhoc(throat);
T_throat(j)=Tc(throat);
pc_throat(j)=pc(throat);
Mc_throat(j)=Mc(throat);
cons_mass_flow_throat(j)=cons_mass_flow(throat);
Vc_throat(j)=Vc(throat);

figure(2)
if j==50
plot(x,cons_mass_flow,'b','linewidth',1.5);
hold on
elseif j==100
plot(x,cons_mass_flow,'r','linewidth',1.5)
hold on
elseif j==150
plot(x,cons_mass_flow,'c','linewidth',1.5)
hold on
elseif j==300
plot(x,cons_mass_flow,'g','linewidth',1.5)
hold on
elseif j==600
plot(x,cons_mass_flow,'y','linewidth',1.5)
hold on
elseif j==900
plot(x,cons_mass_flow,'*','linewidth',1.5)
hold on
elseif j==1400
plot(x,cons_mass_flow,'m','linewidth',1.5)
hold on


end
end
xlabel('Non dimensioanl distance through nozzle(x)','fontsize',10)
ylabel('Non dimensional mass flow rate through nozzle x','fontsize',10);
title('variation of mass flow rate (conservation form) along x at different time steps');
legend({'50 time steps','100 time steps','150 time steps','300 time steps','600 time steps','900 time steps','1400 time steps'});
end

Non conservative :

Conservative :

if the Number of grid points is found to be 61 then,

n=61

non conservative

Grid Dependence :

This test is to determine all the variable ratios at different grid points mainly at n=31 and n=61

CFL effect on non - conservative form for 31 grid points:

The solution blew at CFL value equal to 1.18

CFL data effect conservative form for 31 grid points:

Results :

The CFL value does not have any effect on conservative form primate variables ratio.

From table it can be seen that as we increase the number of grid points we could say that the value approaches very close to the analytical solution.But however the increase in number of grid points increase the simulation time required to initiate the program.

• The time taken for simulation in case of conservative form is low when compared to non conservative form as shown in the simulation table .

• In grid dependence test we observe that the accuracy in case of non consevative form is more as comapared to conservative form as it is very close to exact analytical solution.

• The conservative form provides more accuracy incase of mass flow rate whereas non conservative provide more accuracy in case of primitive variables .However in general we can say that non conservative form has higher accuracy than conservative form.

• No For the first table of cfl effects signifies about the effect on non conservative form.It is observed that as CFL value increases the density,pressure and Temperature ratio increases.But however mach number decreases.

• For the cfl value 1.2 the sloution gets blown off . No wonder it is the threshold CFL value .

• In the second table it represents data for conservative form the increase or decrease in cfl value doesnt effect the solution itself