The objective of this project is to write code to solve the 1D supersonic nozzle flow equations using the MacCormack Method for CFD analysis
Title of Project:
The objective of this project is to write code to solve the 1D supersonic nozzle flow equations using the MacCormack Method for CFD analysis.
- we will be implementing both the conservative and non-conservative forms of the governing equations
we will also perform a grid dependence test by varying grid points along the length.
Introduction:
The nozzle is nothing but a geometry with convergent, divergent and throat section at joint of convergent and divergent portions.
It is an arrangement for converting pressure energy of fluid into kinetic energy. As we know that, with the convergent section, we can increase speed up to a speed of sound and that is known as subsonic zone. whereas with this convergent-divergent nozzle, we can achieve a speed of fluid above the speed of sound and that zone would be supersonic zone. This arrangement is useful for rocket exhaust and certain fluid flow requirements in the industry.
This physical phenomenon can be well described with Navier stokes partial differential equation set of continuity, momentum and energy equations. In CFD analysis these equations are numerically discretized by finite difference approach and further represented in either conservation form(where variables are inside derivative sign), where control volume stood still and fluid moves through it and Non-conservation form(where variables undergo chain rule and separated), where control volume moves with fluid tracking same fluid particle over the time span of observation. then these equations are solved by the MacCormack method, in which the first primitive variables change over time is mapped with forwarding difference approach over space variable and then these updated values change over time is further mapped in backward differencing approach for space variables. This is then followed by averaging both derivatives term is a predictor stage and a corrector stage and then final solution update is done. In this way, this loop continues until we achieve steady-state, where mass flow rate doesn\'t change much over the length of the nozzle and remains within a specified tolerance zone.
Here I have solved these equations by both Non-conservative form solver and conservative form solver, so equations for both forms can be found in chapter 7, John D Anderson jr book named, \'computational fluid dynamics: the basic and applications\'.
Script for MacCormack method calling both Non-conservative form solver and conservative form solver as follows:
clc;
clear all;
close all;
%take input for number of grid points along x axis
n=31;
x=linspace(0,3,n);
dx=x(2)-x(1);
gamma=1.4;
%here we take Courant fredriech lewy number to be 0.5, which is supposed to
%be lesser than 1
CFL=0.5;
%Area profile
A=1+2.2*(x-1.5).^2;
%lets calculate array value number at which A/A*=1,that is throat
n_throat=find(A==1);
error1=9e9;
error2=9e9;
tolerance=1e-5;
%here lets solve 1_D supersonic nozzle flow using macormack method
%for non conservative form solver
tic
[rho,v,T,p,M,mass_flow_rate,rho_throat,T_throat,p_throat,M_throat,k]=non_conservative_solver(n,x,dx,gamma,A,n_throat,CFL,error1,tolerance);
end_time=toc;
fprintf(\'minimum number of iterations for non conservative form:%f\',k);
fprintf(\'time for execution by non conservative solver:%f\',end_time);
%steady profiles
figure(2)
subplot(4,1,1);
plot(x,rho,\'m\');
title(\'density variation along length with time non con\');
legend(\'density ratio\');
grid minor
subplot(4,1,2);
plot(x,p,\'b\');
title(\'pressure variation along length with time non con\');
legend(\'pressure\');
grid minor
subplot(4,1,3);
plot(x,T,\'g\');
title(\'temperature variation along length with time non con\');
legend(\'Temperature\');
grid minor
subplot(4,1,4);
plot(x,M,\'r\');
title(\'mach number variation at throat with time non con\');
legend(\'mach number\');
grid minor
%conservative form solver
tic
[rho_con,v_con,T_con,p_con,M_con,mass_flow_rate_con,rho_throat_con,T_throat_con,p_throat_con,M_throat_con,k1]=conservative_solver(n,x,dx,gamma,A,n_throat,CFL,error2,tolerance);
end_time1=toc;
fprintf(\'minimum number of iterations for conservative form:%f\',k1);
fprintf(\'time for execution by conservative solver:%f\',end_time1);
%steady profiles
figure(4)
subplot(4,1,1);
plot(x,rho_con,\'m\');
title(\'density variation along length with time con\');
legend(\'density ratio\');
grid minor
subplot(4,1,2);
plot(x,p_con,\'b\');
title(\'pressure variation along length with time con\');
legend(\'pressure\');
grid minor
subplot(4,1,3);
plot(x,T_con,\'g\');
title(\'temperature variation along length with time con\');
legend(\'Temperature\');
grid minor
subplot(4,1,4);
plot(x,M_con,\'r\');
title(\'mach number variation at throat with time con\');
legend(\'mach number\');
grid minor
%steady mass flow rate profile for non conservation form and conservation
%form
figure(5)
plot(x,mass_flow_rate)
hold on
plot(x,mass_flow_rate_con)
title(\'steady mass flow rate\')
legend(\'Non conservative\',\'conservative\')
grid minorNon_Conservative form solver:
function [rho,v,T,p,M,mass_flow_rate,rho_throat,T_throat,p_throat,M_throat,k]=non_conservative_solver(n,x,dx,gamma,A,n_throat,CFL,error1,tolerance)
%initial values of rho(density), temperature and velocity are as follows
rho=1-0.3146*x;
T=1-0.2314*x;
v=(0.1+1.09*x).*T.^0.5;
%lets start simulation from t=0 and it will be through while loop which
%checks error1 then proceed for next time step as it will signifies steady
%state of mass flow rate.
k=0;
%let\'s calculate velocity of sound v_sound=(Temperature)^0.5
%let\'s store one copy of values as old
%tfinal=1399*dt;
%t=linspace(0,dt,tfinal);
while(error1>tolerance)
k=k+1;
rho_old=rho;
T_old=T;
v_old=v;
v_sound=(T_old).^0.5;
mass_flow_rate_old=rho_old.*v_old.*A;
Maximum_velocity=v_old+v_sound;
dt_all=CFL*dx./(Maximum_velocity);
dt=min(dt_all);
%predictor time step
for i=2:n-1
dvdx_p=(v(i+1)-v(i))/dx;
dlnAdx_p=(log(A(i+1))-log(A(i)))/dx;
drhodx_p=(rho(i+1)-rho(i))/dx;
dTdx_p=(T(i+1)-T(i))/dx;
%lets assemble in continuity equation
drhodt_p(i)=-rho(i)*dvdx_p-rho(i)*v(i)*dlnAdx_p-v(i)*drhodx_p;
%momentum equation
dvdt_p(i)=-v(i)*dvdx_p-(1/gamma)*(dTdx_p+T(i)/rho(i)*drhodx_p);
%energy equation
dTdt_p(i)=-v(i)*dTdx_p-(gamma-1)*T(i)*(dvdx_p+v(i)*dlnAdx_p);
%Solution update
rho(i)=rho(i)+drhodt_p(i)*dt;
v(i)=v(i)+dvdt_p(i)*dt;
T(i)=T(i)+dTdt_p(i)*dt;
end
%corrector step
for j=2:n-1
dvdx_c=(v(j)-v(j-1))/dx;
dlnAdx_c=(log(A(j))-log(A(j-1)))/dx;
drhodx_c=(rho(j)-rho(j-1))/dx;
dTdx_c=(T(j)-T(j-1))/dx;
%lets assemble in continuity equation
drhodt_c(j)=-rho(j)*dvdx_c-rho(j)*v(j)*dlnAdx_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)*dlnAdx_c);
end
%lets take average of derivative obtained in predictor and corrector
%stage by macormack method
drhodt_avg=(drhodt_p+drhodt_c)*0.5;
dvdt_avg=(dvdt_p+dvdt_c)*0.5;
dTdt_avg=(dTdt_p+dTdt_c)*0.5;
for m=2:n-1
rho(m)=rho_old(m)+drhodt_avg(m)*dt;
v(m)=v_old(m)+dvdt_avg(m)*dt;
T(m)=T_old(m)+dTdt_avg(m)*dt;
end
%lets define boundary conditions
%inlet
rho(1)=1;
v(1)=2*v(2)-v(3);
T(1)=1;
%outlet
rho(n)=2*rho(n-1)-rho(n-2);
v(n)=2*v(n-1)-v(n-2);
T(n)=2*T(n-1)-T(n-2);
%lets calxulate mach number and mass flow rate of flow
M=v./sqrt(T);
mass_flow_rate=rho.*v.*A;
p=rho.*T;
%transient throat values at throat section of nozzle
rho_throat(k)=rho(n_throat);
M_throat(k)=M(n_throat);
T_throat(k)=T(n_throat);
p_throat(k)=p(n_throat);
%lets check convergence
error1=max(abs((mass_flow_rate-mass_flow_rate_old)));
%lets plot data for transient flow
figure(1);
hold on
if k==1
plot(x,mass_flow_rate)
elseif k==10
plot(x,mass_flow_rate)
elseif k==100
plot(x,mass_flow_rate)
elseif k==250
plot(x,mass_flow_rate)
elseif k==500
plot(x,mass_flow_rate)
elseif k==700
plot(x,mass_flow_rate)
end
title(\'transient mass flow rate for non conservation form\');
xlabel(\'L\');
ylabel(\'mass flow rate\');
legend(\'t=1*dt\',\'50*dt\',\'100*dt\',\'200*dt\',\'500*dt\',\'700*dt\');
axis([0 3 0 1.8]);
grid minor
end
Conservative form solver:
function [rho_con,v_con,T_con,p_con,M_con,mass_flow_rate_con,rho_throat_con,T_throat_con,p_throat_con,M_throat_con,k1]=conservative_solver(n,x,dx,gamma,A,n_throat,CFL,error2,tolerance)
%initial values of rho(density), temperature and velocity are as follows
for i=1:length(x)
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);
elseif x(i)>=1.5 && x(i)<=3
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);
p_con=rho_con.*T_con;
%lets define initial solution vectors from primitive variables
U1=rho_con.*A;
U2=rho_con.*v_con.*A;
U3=rho_con.*A.*((T_con./(gamma-1))+(gamma/2)*v_con.^2);
%lets start simulation from t=0 and it will be through while loop which
%checks error1 then proceed for next time step as it will signifies steady
%state of mass flow rate.
k1=0;
while error2>tolerance
k1=k1+1;
mass_flow_rate_old_con=rho_con.*v_con.*A;
U1_old=U1;
U2_old=U2;
U3_old=U3;
%lets define flux vector as function of solution vector and source
%terms as combination of primitive term and solution vector
F1=U2;
F2=(U2.^2./U1)+((gamma-1)/gamma)*(U3-(gamma*U2.^2)./(2*U1));
F3=(gamma*U2.*U3./U1)-(gamma*(gamma-1)/2)*(U2.^3./U1.^2);
%let\'s calculate velocity of sound v_sound=(Temperature)^0.5
v_sound=(T_con).^0.5;
%set time step based on CFL number
Maximum_velocity=v_con+v_sound;
dt_all=CFL*dx./(Maximum_velocity);
dt=min(dt_all);
%predictor time step
for j=2:n-1
%defining source term
J2(j)=(1/gamma).*rho_con(j).*T_con(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(j);
%energy equation
dU3dt_p(j)=-(F3(j+1)-F3(j))/dx;
%update predictor step solution
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
%here value of updated primitive variable and flux vectors is as
%follows
rho_con=U1./A;
v_con=U2./U1;
T_con=(gamma-1)*(U3./U1-(gamma/2)*v_con.^2);
p_con=rho_con.*T_con;
F1=U2;
F2=(U2.^2./U1)+((gamma-1)/gamma)*(U3-(gamma*U2.^2)./(2*U1));
F3=(gamma*U2.*U3./U1)-(gamma*(gamma-1)/2)*(U2.^3./U1.^2);
%corrector derivative step
for j=2:n-1
%defining source term
J2(j)=(1/gamma).*rho_con(j).*T_con(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(j);
%energy equation
dU3dt_c(j)=-(F3(j)-F3(j-1))/dx;
end
%average derivatives are as follows
dU1dt_avg=0.5*(dU1dt_p+dU1dt_c);
dU2dt_avg=0.5*(dU2dt_p+dU2dt_c);
dU3dt_avg=0.5*(dU3dt_p+dU3dt_c);
%final solution stage
for j=2:n-1
U1(j)=U1_old(j)+dU1dt_avg(j)*dt;
U2(j)=U2_old(j)+dU2dt_avg(j)*dt;
U3(j)=U3_old(j)+dU3dt_avg(j)*dt;
end
%lets set boundary conditions
%inlet
rho_con(1)=1;
T_con(1)=1;
U1(1)=rho_con(1)*A(1);
U2(1)=2*U2(2)-U2(3);
U3(1)=U1(1)*(T_con(1)/(gamma-1)+(gamma/2)*v_con(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);
%calculation of primitive variable
rho_con=U1./A;
v_con=U2./U1;
T_con=(gamma-1)*(U3./U1-(gamma/2)*v_con.^2);
p_con=rho_con.*T_con;
%calculation of Mach number, mass flow rate and pressure for
%conservative flow
M_con=v_con./sqrt(T_con);
mass_flow_rate_con=rho_con.*v_con.*A;
p_con=rho_con.*T_con;
% transient throat values
rho_throat_con(k1)=rho_con(n_throat);
T_throat_con(k1)=T_con(n_throat);
p_throat_con(k1)=p_con(n_throat);
M_throat_con(k1)=M_con(n_throat);
%convergence assurance
error2=max(abs(mass_flow_rate_con-mass_flow_rate_old_con));
%transition flow rate profile
figure(3)
hold on
if k1==1
plot(x,mass_flow_rate_con)
elseif k1==50
plot(x,mass_flow_rate_con)
elseif k1==100
plot(x,mass_flow_rate_con)
elseif k1==200
plot(x,mass_flow_rate_con)
elseif k1==500
plot(x,mass_flow_rate_con)
elseif k1==700
plot(x,mass_flow_rate_con)
end
title(\'transient mass flow rate for conservation form\');
xlabel(\'L\');
ylabel(\'mass flow rate conservative flow\');
legend(\'t=1*dt\',\'50*dt\',\'100*dt\',\'200*dt\',\'500*dt\',\'700*dt\');
axis([0 3 0.5 0.8]);
grid minor
end
Solution Results:
- Let\'s figure out the minimum number of cycles for which the simulation should be run in order for convergence and here is output from code for a number of iteration and time for execution as follows.
- The minimum number of iterations for non-conservative form:767.000000
- Time for execution by non-conservative solver:131.140124 Sec
- The minimum number of iterations for conservative form:951.000000
- Time for execution by conservative solver:174.765540 Sec
Results For Non-Conservative form solver for transient mass flow rate along the length of the nozzle and steady-state variation of the primitive variable along the length as follows:
figure(1)
figure(2)
Results For Conservative form solver for transient mass flow rate along the length of the nozzle and steady-state variation of the primitive variable along the length as follows:
figure(3)
figure(4)
Here are graphs showing normalized mass flow rate of fluid across the length of the nozzle for both conservative form and non-conservative form of solver and it is as follows:
figure(5)
Conclusions:
1) Here based on the number of iterations, it looks like a nonconservative solver requires a minimum number of iterations to achieve a steady-state than that of conservative form solver and thus non-conservative form solver is faster than conservative form solver.
2) Here from figure 1, it can be seen that though non-conservative form solver reaches a steady-state faster, this transition is coarser than that of conservative form solver as seen from figure 3.
3) As seen in figure 2 and figure 4, we can say that pressure, temperature and density vary over the length of the nozzle as expected decreases whereas velocity is on the rise as seen from Mach number which is nothing but the ratio of fluid velocity to sound velocity.
4) From figure 5, we have validated that non-conservative form solver is faster but simultaneously it is coarser than conservative form solver.
5) Grid dependence Test:
here we have calculated solutions for varying grid points along the length of the nozzle and we have two cases
case 1: n=31
average mass flow rate non-conservative form solver:
| 0.587925 |
average mass flow rate conservative form solver:
| 0.585173 |
case 2:n=61
average mass flow rate non-conservative form solver:
| 0.583698 |
average mass flow rate conservative form solver:
| 0.58285 |
Here we observed all variables at the throat section of the nozzle as at which value of Mach number should be 1 and thus we compare all of them at this benchmark.
Density ratio Temperature ratio pressure ratio Mach number
Analy- 0.634 0.833 0.528 1
tical
Non-conser- 0.638 0.8364 0.5341 0.9993
vative form
n=31
conser- 0.650 0.8401 0.5466 0.9813
vative form
n=31
Non-conser- 0.6377 0.8355 0.5328 0.9998
vative form
n=61
conser- 0.6385 0.8353 0.5333 0.9992
vative form
n=61
It seems like the non-conservative form is more accurate even with 31 grid points than a conservative form with 61 grid points. Thus Non-conservative form is a grid size independent for given accuracy for mass flow rate.
Final Conclusion:
1) the conservative form is better at mass flow distribution over the length of the nozzle and thus gives better mass conservation.
2) conservative form is much slower than Non-conservative form due to computation effort for calculating primitive variables from flux variables.
3) The non-conservative form has much lesser residual and thus accuracy is being achieved faster than conservative form.