1D SUBSONIC-SUPERSONIC ISENTROPIC NOZZLE FLOW SIMULATION USING MACCORMACK METHOD
PROJECT OBJECTIVE: To study the aspects of a steady isentropic flow analysis through a convergent-divergent nozzle.
ABSTRACT:
In this particular project, a time marching finite-difference solution is set up for the quasi-1D-isentropic nozzle flow simulation.The simulation is being performed by constructing a simplified engineering model of the nozzle.
It is assumed that the flow at the inlet to the nozzle comes from a reservoir .The cross-section of the reservoir is considered to be large and hence the velocity is very small at the inlet.The flow expands isentropically to supersonic speeds at the nozzle exit.The flow is locally subsonic at the convergent section of the nozzle, sonic at the throat and supersonic at the divergent section.
Here, an assumption is made that the flow properties vary only with x-direction, hence termed as quasi 1D flow.
The simulation is performed using both conservation and non-conservation forms of the governing equations.
In real life, there isn't actually much difference between conservation and non-conservation equation, but mathematically,in the non-conservation eqn, we imagine a control volume inside the fluid flow moving with the flow and the parameters of the flow inside the control volume changing in space accordingly.
In conservation form of governing eqn, we imagine the control volume to be fixed in a particular location inside the flow and the flow is passing through it and the parameter of the flow are changing with respect to space.
The derived equations used in this project can be found in computational fluid dynamics book like CFD BY J.D.ANDERSON.
PROGRAMING IN MATLAB:
The non-conservation part and conservation part are initially programmed separately to reduce the complexity in viewing and understanding the code and the flow parameters are computed separately.
While comparing some flow parameters between the non-conservation set and the conservation set, functions are being used to call the parameters to ease the viewing of the results.
An important numerical numerical distinction between the non-conservation and conservation form of the equations is realted to the problems dealing with shock-capturing.
NON-CONSERVATION FORM:
PROGRAM SNIPPET FOR THE NON-CONSERVATION FORM:
% McCormack Method for solving 1D Supersonic Nozzle Flow
% Non-consevation code
clear all
close all
clc
%inputs
n=31;
x=linspace(0,3,n);
dx=x(2)-x(1);
gamma=1.4;
%Data for CFL time loop and convergence
c=0.5;
%calculate initial profile
rho=1-0.3146*x;
t=1-0.2314*x;
v=(0.1+1.09*x).*t.^0.5;
%area profile
a=1+2.2*(x-1.5).^2;
%time step
nt=1000;
%outer time loop
for k=1:nt
rho_old =rho; %retaining copy of the solver variable
v_old=v;
t_old=t;
a_old=sqrt(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;
% dt calculation from CFL number
dt=min(c*dx./(a_old+v_old));
%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);
%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 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
%compute the average time derivative
drhodt=0.5*(drhodt_p+drhodt_c);
dvdt=0.5*(dvdt_p+dvdt_c);
dtdt=0.5*(dtdt_p+dtdt_c);
%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
%applying the boundary condition
%inlet
rho(1)=1;
v(1)=2*v(2)-v(3);
t(1)=1;
%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 variables
M=v./sqrt(t); % Mach number
massflow=rho.*v.*a; % Non dimensional mass flow rate
p=rho.*t; % non dimensional pressure
%Throat parameters
rho_throat(k)=rho(16); %density at throat
t_throat(k)=t(16); %temperature at throat
p_throat(k)=p(16); %pressure at throat
M_throat(k)=M(16); %Mach number at throat
% TO CHECK THE MASS FLOW RATE PROFILE WITH RESPECT TO TIME
figure(1) % mass flow rate plot
hold on
if k==1
plot(x,massflow,'Color','k')
elseif k==50
plot(x,massflow,'color','g')
elseif k==100
plot(x,massflow,'Color','m')
elseif k==150
plot(x,massflow,'color','b')
elseif k==200
plot(x,massflow,'Color','c')
elseif k==700
plot(x,massflow,'Color','y','LineWidth',1.4)
plot(x(end),massflow(end),'color','k','Marker','o')
end
title('Mass Flow rate- Transient State');
xlabel('Non-dimensional Mass Flow')
ylabel('Non-dimensional distance through nozzle')
legend('0\Deltat','50\Deltat','100\Deltat','150\Deltat','200\Deltat','700\Deltat')
axis([0 3 0 2]);
% FLOW PARAMETERS AT THE NOZZLE THROAT
figure(2)
subplot(4,1,1)
plot(rho_throat,'color','b')
title('Time wise variation of flowparameters at the nozzle throat');
ylabel('rho/rho_{o}')
legend('Density ratio')
axis([0 k 0.5 1]);
grid minor;
subplot(4,1,2)
plot(t_throat,'color','m')
ylabel('T/T_{o}')
legend('Temperature ratio')
axis([0 k 0.5 1]);
grid minor;
subplot(4,1,3)
plot(p_throat,'color','g')
ylabel('p/p_{o}')
legend('pressure ratio')
axis([0 k 0.5 1]);
grid minor;
subplot(4,1,4)
plot(M_throat,'color','c')
ylabel('M')
legend('Mach number')
axis([0 k 0.9 1.3]);
grid minor;
xlabel('Number of timesteps')
end
% QUALITATIVE PARAMETERS OF 1 DIMENSIONAL NOZZLE FLOW
figure(3)
subplot(4,1,1)
plot(x,M,'color','c')
ylabel('M')
legend('Mach number','location','northwest')
title('Qualitaive aspects of quasi 1D nozzle flow');
grid minor;
subplot(4,1,2)
plot(x,p,'color','b')
ylabel('p/p_{o}')
legend('pressure ratio','location','southwest')
grid minor;
subplot(4,1,3)
plot(x,rho,'color','r')
ylabel('rho/rho_{o}')
legend('density ratio','location','southwest')
grid minor;
subplot(4,1,4)
plot(x,t,'color','c')
ylabel('T/T_{o}')
legend('Temperature ratio','location','southwest')
grid minor;
TIMEWISE VARIATIONS AT THE THROAT:
The variations of the properties at the throat with the change of iterations are shown in the following plots.It can be seen that the stability in flow simulation has been achieved after nearly 600 iterations.
T
It is observed that the nondimensional density at the throat has changed from 0.528 to 0.531, a 0.57% change over one time step.However in the approach towards the steady state solution, at larger values of time, the changes in the flow field variables from one time step to the next become smaller and zero in the limit of large time.
The largest changes are seen to take place at early times after which final steady values are approached almost asymptotically.It can also be mentioned as the "rubber- band effect".
At early times, the oscillations are large in values and are because these oscillations are associated with various unsteady compression and expansion waves which propagate through the nozzle during the transient process.
MASS FLOW VARIATION WITH ITERATIONS:
The nondimensional mass flow is plotted as function of nondimensioanl distance through the nozzle.Six different curves are shown,each for a different time during the course of the time marching procedure. After 200 time steps, the mass flow distribution is beginning to settle down and after 700 time steps the mass flow distribution is a straight horizental line across the graph, which says that the mass flow value has converged to a constant, steady-state value throughout the nozzle.
TIME FOR SIMULATION:
STEADY STATE PROFILE:
After the simulation is over, the staedy state is reached and the steady state profile for pressure ratio, density ratio, temperature ratio and mach number are plotted along the domain of the nozzle.
The results seen are almost identical to the analytical solution and as seen the density, temperature and pressure ratio are dropping as the fluid is accelerating along the domain.As mention earlier in the abstract that the flow is initially subsonic will reach sonic stateat the nozzle throat and will become supersonic at the exit, it is eventually verified in the plot of the mach number along the domain of the nozzle.
CONSERVATION FORM:
PROGRAM SNIPPET FOR THE CONSERVATION FORM:
clear all
close all
clc
%inputs
n=31;
x=linspace(0,3,n);
dx=x(2)-x(1);
gamma=1.4;
c=0.5;
%area
a=1+2.2*(x-1.5).^2;
for i=1:n
%calculate initial profile
if(x(i)>=0 && 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);
else
rho(i)=0.634-0.3879*(x(i)-1.5);
t(i)=0.833-0.3507*(x(i)-1.5);
end
end
v=0.59./(rho.*a);
%calculation of generic elements of the solution vector
u1=rho.*a;
u2=rho.*a.*v;
u3=rho.*a.*(((1/(gamma-1))*t)+((gamma/2)*(v.^2)));
%time step
nt=1400;
%outer time loop
for k=1:nt
u1_old =u1; %retaining copy of the solver variable
u2_old=u2;
u3_old=u3;
%dt calculation from CFL number
dt=min(c*dx./(sqrt(t)+v));
%defining flux vectors
f21=((u2).^2)./(u1);
f31=((u2).*(u3))./(u1);
f32=((u2).^3)./((u1).^2);
%predictor method
for j=2:n-1
df1dx=(u2(j+1)-u2(j))/dx;
df2dx=((f21(j+1)+(1/gamma)*(gamma-1)*(u3(j+1)-(gamma/2)*f21(j+1)))-(f21(j)+(1/gamma)*(gamma-1)*(u3(j)-(gamma/2)*f21(j))))/dx;
df3dx=((gamma*f31(j+1)-0.5*gamma*(gamma-1)*f32(j+1))-(gamma*f31(j)-0.5*gamma*(gamma-1)*f32(j)))/dx;
dlogadx=(log(a(j+1))-log(a(j)))/dx;
J2=(gamma-1)/gamma*(u3-(gamma/2)*f21)*dlogadx;
%continuity equation
du1dt_p(j)=-df1dx;
%momentum equation
du2dt_p(j)=-df2dx+J2(j);
%energy equation
du3dt_p(j)=-df3dx;
%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
rho = u1./a;
v = u2./u1;
t = (gamma-1).*((u3./u1)-(gamma.*(v.^2)/2));
%defining flux vectors
f21=((u2).^2)./(u1);
f31=((u2).*(u3))./(u1);
f32=((u2).^3)./(u1).^2;
%corrector method
for j=2:n-1
df1dx=(u2(j)-u2(j-1))/dx;
df2dx=((f21(j)+(1/gamma)*(gamma-1)*(u3(j)-(gamma/2)*f21(j)))-(f21(j-1)+(1/gamma)*(gamma-1)*(u3(j-1)-(gamma/2)*f21(j-1))))/dx;
df3dx=((gamma*f31(j)-0.5*gamma*(gamma-1)*f32(j))-(gamma*f31(j-1)-0.5*gamma*(gamma-1)*f32(j-1)))/dx;
dlogadx=(log(a(j))-log(a(j-1)))/dx;
J2=(gamma-1)/gamma*(u3-(gamma/2)*f21)*dlogadx;
%continuity equation
du1dt_c(j)=-df1dx;
%momentum equation
du2dt_c(j)=-df2dx+J2(j);
%energy equation
du3dt_c(j)=-df3dx;
end
%compute the average time derivative
du1dt=0.5*(du1dt_p+du1dt_c);
du2dt=0.5*(du2dt_p+du2dt_c);
du3dt=0.5*(du3dt_p+du3dt_c);
%final solution update
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
%applying the boundary condition
%inlet
u2(1)=2*u2(2)-u2(3);
%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);
%calculating the flow parameters
rho=u1./a;
v=u2./u1;
t=(gamma-1)*((u3./u1)-((gamma/2)*v.^2));
%calculation of variables
M=v./sqrt(t); % Mach number
massflow=rho.*v.*a; % Non dimensional mass flow rate
p=rho.*t; % non dimensional pressure
%Throat parameters
rho_throat(k)=rho(16); %density at throat
t_throat(k)=t(16); %temperature at throat
p_throat(k)=p(16); %pressure at throat
M_throat(k)=M(16); %Mach number at throat
% TO CHECK THE MASS FLOW RATE PROFILE IN CONSERVATION FORM WITH RESPECT TO TIME
figure(1) % mass flow rate plot
hold on
if k==1
plot(x,massflow,'Color','k','LineWidth',1.4)
elseif k==100
plot(x,massflow,'Color','m','LineWidth',1.4)
elseif k==200
plot(x,massflow,'Color','c','LineWidth',1.4)
elseif k==700
plot(x,massflow,'Color','y','LineWidth',1.4)
plot(x(end),massflow(end),'color','k','Marker','*')
end
title('Mass Flow rate- Transient State');
xlabel('Non-dimensional Mass Flow')
ylabel('Non-dimensional distance through nozzle')
legend('0\Deltat','100\Deltat','200\Deltat','700\Deltat')
axis([0 3 0.5 0.8]);
% FLOW PARAMETERS AT THE NOZZLE THROAT IN CONSERVATION FORM
figure(2)
subplot(4,1,1)
plot(rho_throat,'color','b','LineWidth',1.4)
title('Time wise variation of flowparameters at the nozzle throat');
ylabel('rho/rho_{o}')
legend('Density ratio')
axis([0 k 0.5 1]);
grid minor;
subplot(4,1,2)
plot(t_throat,'color','m','LineWidth',1)
ylabel('T/T_{o}')
legend('Temperature ratio')
axis([0 k 0.6 1.2]);
grid minor;
subplot(4,1,3)
plot(p_throat,'color','g','LineWidth',1)
ylabel('p/p_{o}')
legend('pressure ratio')
axis([0 k 0.4 1]);
grid minor;
subplot(4,1,4)
plot(M_throat,'color','c','LineWidth',1)
ylabel('M')
legend('Mach number')
axis([0 k 0.8 1.3]);
grid minor;
xlabel('Number of timesteps')
end
% QUALITATIVE PARAMETERS OF CONSERVATION FORM OF 1 DIMENSIONAL NOZZLE FLOW
figure(3)
subplot(4,1,1)
plot(x,M,'color','c','LineWidth',1)
ylabel('M')
legend('Mach number','location','northwest')
title('Qualitaive aspects of quasi 1D nozzle flow');
grid minor;
subplot(4,1,2)
plot(x,p,'color','b','LineWidth',1)
ylabel('p/p_{o}')
legend('pressure ratio','location','southwest')
grid minor;
subplot(4,1,3)
plot(x,rho,'color','r','LineWidth',1)
ylabel('rho/rho_{o}')
legend('density ratio','location','southwest')
grid minor;
subplot(4,1,4)
plot(x,t,'color','c','LineWidth',1)
ylabel('T/T_{o}')
legend('Temperature ratio','location','southwest')
grid minor;
TIMEWISE VARIATION OF FLOW PARAMETERS AT THE THROAT:
The properties seem to be almost stable even in early iterations and is seen to reach stability soon. The convergence is looked for in the plot.
MASSFLOW RATE VARIATIONS:
The nondimensional mass flow is plotted versus x/L for for several different values of time.At 0 time step, the initial conditions are represented and the transient flow is seen to deviate away from the initial conditions.The results after 100 time steps shows a somewhat of a "humped" distribution.After 200 time steps, the mass flow distribution is becoming more constant, and after 700 time steps, it is almost equal to constant value.Moreover, it is quite close to the exact analytical value of 0.579.
COMPARISON OF CONSERVATION AND NONCONSERVATION MASSFLOW AT VARIOUS TIME:
Comparing the corresponding results of conservation forem of equations with that obtained with the nonconserbation form of equations, it is seen that in conservation for, the variations in mass flow are much less severe.It is supposed that the milder behaviour is predominantly due to assumed initial conditions of constant mass flow.
TIME FOR SIMULATION:
STEADY STATE PROFILE:
After the simulation is over, the staedy state is reached and the steady state profile for pressure ratio, density ratio, temperature ratio and mach number are plotted along the domain of the nozzle.
The steady state results obtained from the time marching solution of the governing equations in conservation form are essentially the same as those used in nonconservation form.The results are obtained after 1400 time steps. Hence it can be concluded that in all practical purposes both forms of governing equations give the same results.
COMPARISON OF STEADY MASSFLOW VARIATIONS:
In this comparison, the scale of ordinate in massflow is greatly magnified.It can be seen that the steady state mass flow distribution predicted by the conservation form of equations is much more satisfactory than that obtained from the non-conservation form.
1. The conservation form gives a distribution that is much closer to being a constant.In contrast, the nonconservation results have a sizeable variation with some spurious oscillations at both the inflow and the outflow boundary.
2. The steady state mass flow results obtained with the conservation form are much closer to the analytical solution though having some variations in the beginning.
THE PROGRAM SNIPPET FOR COMPARISON OF MASSFLOW PROFILES USING FUNCTIONS:
THE MAIN CODE:
clear all
close all
clc
n=31;
x=linspace(0,3,n);
dx=x(2)-x(1);
gamma=1.4;
c=0.5;
a=1+2.2*(x-1.5).^2;
%non conservation equation
tic
[rho,t,p,massflow,M,k]=non_conservation(n,x,dx,gamma,a,c);
time_nc=toc;
%conservation eqn
tic
[rho_con,t_con,p_con,massflow_con,M_con,k1]=conservation(n,x,dx,gamma,a,c);
time_c=toc;
figure(1)
plot(x,massflow,'color','r','linewidth',1.2)
hold on
plot(x,massflow_con,'color','b','linewidth',1.2)
legend('mass flow non-conservation','mass flow conservation')
title('steady mass flow profiles')
NONCONSERVATION FUNCTION:
function[rho,t,p,massflow,M,k]=non_conservation(n,x,dx,gamma,a,c)
%calculate initial profile
rho=1-0.3146*x;
t=1-0.2314*x;
v=(0.1+1.09*x).*t.^0.5;
%time step
nt=1000;
for k=1:nt
rho_old =rho; %retaining copy of the solver variable
v_old=v;
t_old=t;
for j=2:n-1 %predictor
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;
dt=min(c*dx./(sqrt(t)+v));
drhodt_p(j)=-rho(j)*dvdx-rho(j)*v(j)*dlogadx-v(j)*drhodx;
dvdt_p(j)=-v(j)*dvdx-(1/gamma)*(dtdx+(t(j)/rho(j))*drhodx);
dtdt_p(j)=-v(j)*dtdx-(gamma-1)*t(j)*(dvdx+v(j)*dlogadx);
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
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;
drhodt_c(j)=-rho(j)*dvdx-rho(j)*v(j)*dlogadx-v(j)*drhodx;
dvdt_c(j)=-v(j)*dvdx-(1/gamma)*(dtdx+(t(j)/rho(j))*drhodx);
dtdt_c(j)=-v(j)*dtdx-(gamma-1)*t(j)*(dvdx+v(j)*dlogadx);
end
drhodt=0.5*(drhodt_p+drhodt_c);
dvdt=0.5*(dvdt_p+dvdt_c);
dtdt=0.5*(dtdt_p+dtdt_c);
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=v./sqrt(t); % Mach number
massflow=rho.*v.*a; % Non dimensional mass flow rate
p=rho.*t; % non dimensional pressure
endCONSERVATION FUNCTION:
function [rho_con,t_con,p_con,massflow_con,M_con,k1]=conservation(n,x,dx,gamma,a,c)
for i=1:n
if(x(i)>=0 && x(i)<=0.5)
rho_con(i)=1;
t_con(i)=1;
elseif (x(i)>0.5 && x(i)<=1.5)
rho_con(i)=1-0.366*(x(i)-0.5);
t_con(i)=1-0.167*(x(i)-0.5);
else
rho_con(i)=0.634-0.3879*(x(i)-1.5);
t_con(i)=0.833-0.3507*(x(i)-1.5);
end
end
v_con=0.59./(rho_con.*a);
u1=rho_con.*a;
u2=rho_con.*a.*v_con;
u3=rho_con.*a.*(((1/(gamma-1))*t_con)+((gamma/2)*(v_con.^2)));
nt=1400;
for k1=1:nt
u1_old =u1; %retaining copy of the solver variable
u2_old=u2;
u3_old=u3;
dt=min(c*dx./(sqrt(t_con)+v_con));
f21=((u2).^2)./(u1);
f31=((u2).*(u3))./(u1);
f32=((u2).^3)./((u1).^2);
for j=2:n-1
df1dx=(u2(j+1)-u2(j))/dx;
df2dx=((f21(j+1)+(1/gamma)*(gamma-1)*(u3(j+1)-(gamma/2)*f21(j+1)))-(f21(j)+(1/gamma)*(gamma-1)*(u3(j)-(gamma/2)*f21(j))))/dx;
df3dx=((gamma*f31(j+1)-0.5*gamma*(gamma-1)*f32(j+1))-(gamma*f31(j)-0.5*gamma*(gamma-1)*f32(j)))/dx;
dlogadx=(log(a(j+1))-log(a(j)))/dx;
J2=(gamma-1)/gamma*(u3-(gamma/2)*f21)*dlogadx;
du1dt_p(j)=-df1dx;
du2dt_p(j)=-df2dx+J2(j);
du3dt_p(j)=-df3dx;
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
rho_con = u1./a;
v_con = u2./u1;
t_con = (gamma-1).*((u3./u1)-(gamma.*(v_con.^2)/2));
f21=((u2).^2)./(u1);
f31=((u2).*(u3))./(u1);
f32=((u2).^3)./(u1).^2;
for j=2:n-1
df1dx=(u2(j)-u2(j-1))/dx;
df2dx=((f21(j)+(1/gamma)*(gamma-1)*(u3(j)-(gamma/2)*f21(j)))-(f21(j-1)+(1/gamma)*(gamma-1)*(u3(j-1)-(gamma/2)*f21(j-1))))/dx;
df3dx=((gamma*f31(j)-0.5*gamma*(gamma-1)*f32(j))-(gamma*f31(j-1)-0.5*gamma*(gamma-1)*f32(j-1)))/dx;
dlogadx=(log(a(j))-log(a(j-1)))/dx;
J2=(gamma-1)/gamma*(u3-(gamma/2)*f21)*dlogadx;
du1dt_c(j)=-df1dx;
du2dt_c(j)=-df2dx+J2(j);
du3dt_c(j)=-df3dx;
end
du1dt=0.5*(du1dt_p+du1dt_c);
du2dt=0.5*(du2dt_p+du2dt_c);
du3dt=0.5*(du3dt_p+du3dt_c);
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
u2(1)=2*u2(2)-u2(3);
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);
rho_con=u1./a;
v_con=u2./u1;
t_con=(gamma-1)*((u3./u1)-((gamma/2)*v_con.^2));
M_con=v_con./sqrt(t_con); % Mach number
massflow_con=rho_con.*v_con.*a; % Non dimensional mass flow rate
p_con=rho_con.*t_con; % non dimensional pressure
endTHE MATTER OF GRID DEPENDENCE:
The matter of grid dependence is a seious consideration in CFD.The value of flow field variables are found to be quite different with different grid points.We must, if possible increase the number of grid points until we reach a solution which is no longer sensitive to the number of points. This stage is known as grid independence.
Here the grid points are varied from 31 to 61 and the same calculations are performed.The results predict the nonconservation solvers to be more accurate.
MACH NUMBER TEMPERATURE DENSITY PRESSURE
RATIO RATIO RATIO
ANALYTICAL 1.0 0.833 0.634 0.528
NONCONSERVATION-31 0.999 0.836 0.639 0.534
NONCONSERVATION-61 0.987 0.839 0.648 0.544
CONSERVATION-31 0.999 0.835 0.637 0.532
CONSERVATION-61 0.989 0.838 0.644 0.540