Week 7 - Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
AIM:
1.To simulate the isentropic flow through a quasi 1D subsonic-supersonic nozzle.
2.To derive the conservation and non-conservation forms of the governing equations and solve them using the MacCormack's technique.
3.To plot the distribution of primitive variables, mass flow rate through the nozzle and time-wise variations of primitive variables at throat for both the conservative and non-conservative forms.
4.To investigate the difference between the two forms of governing equations by comparing their solutions.
DESCRIPTION:
The subsonic-supersonic nozzle of length L consists of two sections namely the convergent and divergent sections separated by a throat region.
The above nozzle can be set up by
1.The governing equations will be expressed in terms of PDE suitable for time marching solution of quasi 1D flow.
2.The finite difference expressions pertaining to MacCormack's technique is to be set up.
3.Calculation of time steps and treatment of initial and boundary conditions will be formulated.
For fluid flow through nozzle, the flow field variables are expressed in non-dimensional form.By doing so, the numerical values of these variables becomes easier to present and compare the results.
FLOW FIELD VARIABLES IN NON-DIMENSIONAL FORM
`T'= T/T^o`
`rho' = rho/rho^o`
`x' = x/L`
`V' = V/aunderset(o)`
Where `aunderset(o)`is the speed of sound
`ao = sqrt(gammaRTo )`
`t' = t/((L)/(ao))`
`A' = A/A^*`
where `A^*` is throat area
GOVERNING EQUATIONS IN NON-CONSERVATIVE FORM (NON-DIMENSIONAL)
The continuity equation is
`(delrho')/(delt')=-rho'.(rhoV')/(rhox')-rho'V'.(rho(ln A'))/(delx')-V'.(delrho')/(rhox')`
The momentum equation is
`(rhoV')/(rhot')=-V'.(rhoV')/(rhox')-1/gamma.((delT')/(delx')+(T')/(rho').(delrho')/(delx'))`
The energy equation is
`(delT')/(delt')=-V'.(delT')/(delx')-(gamma-a).T'.[(delV')/(delx')+V'.(del(ln A'))/(delx')]`
GOVERNING EQUATIONS IN CONSERVATIVE FORM:
For continuity equation is
`(del(rho'A'))/(delt')+(del(rho'A'V'))/(delx')=0`
`(del(rho'A'V'))/(delt')+(del[rho'A'V'^2+(1/gamma)rho'A'])/(delx')=1/gammarho'(delA')/(delx')`
For energy equation is
`(del[rho'((e')/(gamma-1)+gamma/2V'^2)A'])/(delt')+(del[rho'((e')/(gamma-1)+gamma/2V'^2)V'A'+rho'A'V'])/(delx')=0`
writing generic form for above terms as
`U1=rho'A'`
`U2=rho'A'V'`
`U3=rho'((e')/(gamma-1)+gamma/2V'^2)A'`
`F1=rho'A'V'`
`F2=rho'A'V'^2+1/gammarho'A'`
`F3=rho'((e')/(gamma-1)+gamma/2V'^2)V'A'+rho'A'V'`
`J2=1/gammarho'(delA')/(delx')`
writing the above conservation equations for continuity,momentum and energy with substituting above terms we get as
For continuity
`(delU1)/(delt')=-(delF1)/(delx')`
`(delU2)/(delt')=-(delF2)/(delx')+J2`
`(delU3)/(delt')=-(delF3)/(delx')`
here F1=U2 so we get
`F2=(U2^2)/(U1)+(gamma-1)/gamma(U3-gamma/2(U2^2)/(U1))`
`F3=gamma(U2U3)/(U1)-(gamma(gamma-1))/2(U2^3)/(U1^2)`
`J2=(gamma-1)/(gamma)(U3-gamma/2(U2^2)/(U1))(del(ln A'))/(delx')`
MESHING:
The meshing has been taken from reference from computational fluid dynamics J.D Anderson for grid point distribution along the nozzle.
CALCULATION OF TIMESTEPS:
`Deltati=C.(Deltax)/(ai+Vi)`
Local speed of sound `a=sqrtT`
courant number C=0.5
`Deltat=min(Deltat1,Deltat2,....Deltatn)`
DISCRETIZATION BY MACCORMACK'S METHOD:
For predictor step
`(delf)/(deltp)`=Forward difference method for space terms (for derivatives w.r.t x)
For solution update
`f'=fi+(delf)/(deltp).Deltat`
For corrector step
`(delf')/(deltc)`=Backward difference method for space terms (derivative w.r.t x)
`(delf)/(delt(avg))=0.5.((delf)/(deltp)+(delf)/(deltc))`
For final solution update
`f^2=fi+(delf)/(delt(avg)).Deltat`
BOUNDARY CONDITIONS:
The boundary conditions are taken from Computational fluid dynamics by J.D Anderson.
Flow is unsteady and inviscid which is governed by hyperbolic equations for 1D unsteady flow, there exists two real characteristic lines for a point in a plane.
Boundary condition is provided at the domain of dependence
1.If characteristic lines flow from left to right, DIRICHLET B.C is provided which helps to calculate values in domain of influence.Boundary conditions need to be specified.
2.If characteristic lines flow from right to left, NEUMANN B.C is provided which helps to calculate values at the boundary by using values at the inner nodes.Boundary conditions need to be floated.
The appropriate initial and boundary conditions for both non-conservative and conservative forms are taken from the conputational fluid dynamic book by J.D Anderson.
So, The nozzle area `A=1+2.2.(x-1.5)^2` where `0<=x<=3`
Pressure `P=rho.T`
Mach number `M=V/sqrtT`
Mass flow rate = `rho.A.V`
MAIN CODE FOR NOZZLE FLOW
clear all
close all
clc
%inputs
n=31; %number of grid points
x=linspace(0,3,31);
dx=x(2)-x(1);
tic
[t,rho,p,v,machnumber,mass_flow_rate,rho_throat,v_throat,p_throat,t_throat,mach_throat]=non_conservative(n,x,dx);
time1=toc;
tic
[tc,rho_c,pc,vc,machnumber_c,mass_flow_ratec,rho_throatc,v_throatc,p_throatc,t_throatc,mach_throatc,dt]=conservative(n,x,dx);
time2=toc;
figure(1)
subplot(4,1,1)
plot(x,t)
xlabel('Nozzle length')
ylabel('temperature')
subplot(4,1,2)
plot(x,rho)
xlabel('Nozzle length')
ylabel('Density')
subplot(4,1,3)
plot(x,p)
xlabel('Nozzle length')
ylabel('pressure')
subplot(4,1,4)
plot(x,machnumber)
xlabel('Nozzle length')
ylabel('Machnumber')
title('Non-conservative form of distribution of primitive variables along nozzle length')
figure(2)
legend('1st time step','50th time step','100th time steo','200th time step','700th time step','1400th time step');
title('Non-conservative form:variation of mass flow distribution at different time steps');
xlabel('Nozzle length')
ylabel('Mass flow rate')
figure(3)
subplot(4,1,1)
plot(t_throat,'LineWidth',1.5)
xlabel('Time Iterations')
ylabel('temperature at throat')
subplot(4,1,2)
plot(rho_throat,'LineWidth',1.5)
xlabel('Time Iterations')
ylabel('Density at throat')
subplot(4,1,3)
plot(p_throat,'LineWidth',1.5)
xlabel('Time Iterations')
ylabel('Pressure at throat')
subplot(4,1,4)
plot(mach_throat,'LineWidth',1.5)
xlabel('Time Iterations')
ylabel('Mach number at throat')
title('Non-conservative form for time wise variations of primitive variables at throat')
figure(4)
subplot(4,1,1)
plot(x,tc,'LineWidth',1.5)
xlabel('Nozzle length')
ylabel('temperature')
subplot(4,1,2)
plot(x,rho_c,'LineWidth',1.5)
xlabel('Nozzle length')
ylabel('Density')
subplot(4,1,3)
plot(x,pc,'LineWidth',1.5)
xlabel('Nozzle length')
ylabel('Pressure')
subplot(4,1,4)
plot(x,machnumber_c,'LineWidth',1.5)
xlabel('Nozzle length')
ylabel('Mach number')
title('Conservative form for distribution of primitive variables along nozzle length')
figure(5)
legend('1st time step','50th time step','100th time steo','200th time step','700th time step','1400th time step');
title('conservative form:variation of mass flow distribution at different time steps');
xlabel('Nozzle length')
ylabel('Mass flow rate')
figure(6)
subplot(4,1,1)
plot(t_throatc)
xlabel('Time Iterations')
ylabel('temperature')
subplot(4,1,2)
plot(rho_throatc)
xlabel('Time Iterations')
ylabel('Density')
subplot(4,1,3)
plot(p_throatc)
xlabel('Time Iterations')
ylabel('Pressure')
subplot(4,1,4)
plot(mach_throatc)
xlabel('Time Iterations')
ylabel('Mach Number')
title('Conservative form Time-wise variations of primitive variables at throat')
figure(7)
plot(x,mass_flow_rate,'LineWidth',1.5)
plot(x,mass_flow_ratec,'LineWidth',1.5)
legend('Non-conservative form','Conservative form')
xlabel('Nozzle length')
ylabel('Mass flow rate')
title('Comparision of Normalized mass flow rate distribution of both forms')FUNCTION CODE FOR NON_CONSERVATIVE FORM
function [t,rho,p,v,machnumber,mass_flow_rate,rho_throat,v_throat,p_throat,t_throat,mach_throat]=non_conservative(n,x,dx)
%input conditions initially
rho = 1-0.3146*x;
t=1-0.2314*x;
v=(0.1+1.09*x).*t.^0.5;
gamma = 1.4;
%area
a = 1+2.2.*(x-1.5).^2;
%time steps
nt=1400;
dt=min(0.5*dx./(v+sqrt(t)));
%outer time step
for k=1:nt
%assigning old values
rho_old = rho;
v_old=v;
t_old=t;
%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;
%continuity equation
drhodt_p(j)=-rho(j)*dvdx-rho(j)*v(j)*dlogadx-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)*dlogadx));
%updating solution
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 method
for j=2:n-1
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;
%continuity equation
drhodt_c(j)=-rho(j)*dvdx-rho(j)*v(j)*dlogadx-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)*dlogadx));
end
%Average time derivative
drhodt=(drhodt_p+drhodt_c)/2;
dvdt=(dvdt_p+dvdt_c)/2;
dtdt=(dtdt_p+dtdt_c)/2;
%Final solution update
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
%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);
%Calculation of pressure,mach number and mass flow rate
p=rho.*t;
machnumber=v./sqrt(t);
mass_flow_rate=rho.*a.*v;
%Calculation of pressure,mach number and mass flow rate at throat
p_throat(k)=p(16);
v_throat(k)=v(16);
mach_throat(k)=machnumber(16);
t_throat(k)=t(16);
rho_throat(k)=rho(16);
%variation of mass flow rate distribution in nozzle ar different
%timesteps
figure(2)
if k==1
plot(x,mass_flow_rate,'LineWidth',1.5)
hold on
elseif k==50
plot(x,mass_flow_rate,'cyan','LineWidth',1.5)
elseif k==100
plot(x,mass_flow_rate,'r','LineWidth',1.5)
elseif k==200
plot(x,mass_flow_rate,'m','LineWidth',1.5)
elseif k==700
plot(x,mass_flow_rate,'k','LineWidth',1.5)
elseif k==1400
plot(x,mass_flow_rate,'g','LineWidth',1.5)
end
end
FUNCTION CODE FOR CONSERVATIVE FORM
function [tc,rho_c,pc,vc,machnumber_c,mass_flow_ratec,rho_throatc,v_throatc,p_throatc,t_throatc,mach_throatc,dt]=conservative(n,x,dx)
%Initital profiles
for i=1:n
if(0<=x(i) && x(i)<=0.5)
rho_c(i)=1;
tc(i)=1;
elseif (0.5<=x(i) && x(i)<=1.5)
rho_c(i)=1-0.366*(x(i)-0.5);
tc(i)=1-0.167*(x(i)-0.5);
elseif (1.5<=x(i) && x(i)<=3.5)
rho_c(i)=0.634-0.3879*(x(i)-1.5);
tc(i)=0.833-0.3507*(x(i)-1.5);
end
end
%area
a=1+2.2*(x-1.5).^2;
gamma=1.4;
vc=0.59./(rho_c.*a);
%solution vectors
u1=rho_c.*a;
u2=rho_c.*a.*vc;
u3=rho_c.*(tc/(gamma-1)+(gamma/2)*vc.^2).*a;
%time steps
nt=1400
dt=min(0.5*dx./(vc+sqrt(tc)));
%outer time step
for k=1:nt
%Assigning old values
u1_old=u1;
u2_old=u2;
u3_old=u3;
%flux vectors
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 j=2:n-1
j2=(1/gamma)*rho_c(j)*tc(j)*((a(j+1)-a(j))/dx);
%continuity equation
du1dt_p(j)=-(f1(j+1)-f1(j))/dx;
%momentum equation
du2dt_p(j)=(-(f2(j+1)-f2(j))/dx)+j2;
%energy equation
du3dt_p(j)=-(f3(j+1)-f3(j))/dx;
%solution update
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
%updating primitive variables
rho_c=u1./a;
tc=(gamma-1)*((u3)./(u1)-(gamma/2)*((u2)./(u1)).^2);
%updating flux vectors
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 j=2:n-1
j2=(1/gamma)*rho_c(j)*tc(j)*((a(j)-a(j-1))/dx);
%continuity equation
du1dt_c(j)=-(f1(j)-f1(j-1))/dx;
%momentum equation
du2dt_c(j)=(-(f2(j)-f2(j-1))/dx)+j2;
%energy equation
du3dt_c(j)=-(f3(j)-f3(j-1))/dx;
end
%Average time derivative
du1dt=0.5*(du1dt_c+du1dt_p);
du2dt=0.5*(du2dt_p+du2dt_c);
du3dt=0.5*(du3dt_p+du3dt_c);
%Final solution update
for m=2:n-1
u1(m)=u1_old(m)+du1dt(m)*dt;
u2(m)=u2_old(m)+du2dt(m)*dt;
u3(m)=u3_old(m)+du3dt(m)*dt;
end
%Boundary conditions
%Inlet
u1(1)=a(1);
u2(1)=2*u2(2)-u2(3);
u3(1)=u1(1)*(tc(1))/(gamma-1)+(gamma/2)*(vc(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);
%Final update of primitive variables
rho_c=u1./a;
tc=(gamma-1)*((u3)./(u1)-(gamma/2)*((u2)./(u1)).^2);
vc=(u2)./(u1);
%Calculation of pressure,mach number and mass flow rate
pc=rho_c.*tc;
machnumber_c=vc./sqrt(tc);
mass_flow_ratec=rho_c.*a.*vc;
%Calculation of pressure,mach number and mass flow rate at throat
p_throatc(k)=pc(16);
v_throatc(k)=vc(16);
mach_throatc(k)=machnumber_c(16);
t_throatc(k)=tc(16);
rho_throatc(k)=rho_c(16);
%variation of mass flow rate distribution in nozzle ar different
%timesteps
figure(5)
if k==1
plot(x,mass_flow_ratec,'LineWidth',1.5)
hold on
elseif k==50
plot(x,mass_flow_ratec,'cyan','LineWidth',1.5)
elseif k==100
plot(x,mass_flow_ratec,'r','LineWidth',1.5)
elseif k==200
plot(x,mass_flow_ratec,'m','LineWidth',1.5)
elseif k==700
plot(x,mass_flow_ratec,'k','LineWidth',1.5)
elseif k==1400
plot(x,mass_flow_ratec,'g','LineWidth',1.5)
end
endEXPLAINATION OF CODE:
1.Firstly the grid size n,x and dx are defined.
2.Then the function for non-conservation form is defined with initial profiles (rho,t,v),area,time step.
3.Loop for outer time setp is initiated and Mccormack's method is carried out.
4.Predictor method followed by solution update is done.Corrector method followed by final solution update is done.
5.The Boundary conditions for inlet and outlet are defined and then pressure,mach number and mass flow rates are calculated.
6.The values for primitive variables at the throat are assigned to variables with size k so to run for each time step the throat values are calculated.
7.Mass flow rates for different time steps are plotted.
8.The above process is repeated for 1400 time steps and function returns 1x1400 values for primitive variables to main function.
9.Then the conservation function is executed.The same process followed as non-conservation form except that fluc vectors(F1,F2,F3) and solution vectors (U1,U2,U3) are used.
10.The McCormack's method, the flux vectors and solution vectors are used and them primitive variables are updated.
11.The above process is repeated for 1400 time steps and function returns 1x1400 values for primitive variables to main function.
12.The simulation time for conservation and non-conservation form is recorded using tic-toc for time1, and time2.
13.Then plots are plotted.
RESULTS:
For Non-conservative form,the temperature, density and pressure values decrease as the nozzle length increases but mach number increases as nozzle length increases.
For Non-conservative form,the mass flow rate at time steps 700th and 1400th does not change as nozzle length increases.So we can conclude that it has reached steady state at 700th timestep.
For non-conservative for,there are some fluctuations in primitive variables at the initial time iterations and then the plots become linear which implies that steady state is attained.
For conservative form, the temperature, density and pressure values decrease as the nozzle length increases but mach number increases as nozzle length increases.
For conservative form,the mass flow rates at 700th and 1400th time steps does not change as nozzle length increases.So we can conclude that it has reached steady state at 700th time step.
For conservative form, there are some fluctuations in primitive variables at the initial time iterations and then the plots become linear which implies that steady state is attained.
It can be seen that conservation form yields the better mass flow distribution than non-conservative form.
GRID TEST
When changed the grid size from n=31 to 2n,3n... the primitive variable values are almost same as we got for n=31 for both the forms.So, we conclude that grid size is independent for both the forms.
CONCLUSION:
Thus we can conclude that the isentropic flow through a quasi 1D subsonic-supersonic nozzle is simulated successfully by using governing equations using Maccormack's technique and following observations are noticed
1.the conservation form does better in mass conserving than non-conservation form.
2.conservation form requires extra effort as we need to calculate flux variables and vectors.
3.Simulation time in conservation form is less than non-conservation.