Cfd analysis of 1d supersonic nozzle flow using MacCormack method
Goal: solve the 1d supersonic nozzle flow using Maccormack method
objective:
1) Implementing both conservative and non conservative form of governing equations by writing separate functions.
2) figure out minimum number of cycles for which simulation should be run for convergence.
3) compare normalized mass flow rate from both the conservative and non conservative form.
4) determine which method is faster.
#Main code
clear all
close all
clc
%inputs
n=31;
x=linspace(0,3,n);
xr=x./3; % normalized length
dx=x(2)-x(1);
gamma=1.4;
% area
a=1+2.2*(x-1.5).^2;
%location of throat
n_throat=find(a==1);
%paramter which keeps count of no of cycles for convergence
g=0;
tol=1e-6;
error=9e9;
%Non Conservative Form Solver
tic
[rho,v,t,g,m,p,mach,Mach_throat,Pressure_ratio_throat,Density_ratio_throat,Temperature_ratio_throat,k]=non_conservative(n,x,dx,a,gamma,g,tol,error,n_throat);
t_non_consv=toc;
fprintf(\'minimum number of iterations for convergence of non conservative form:%f\',g);
fprintf(\'time taken by non conservstive solver:%f\',t_non_consv);
%steady state profiles for non conservative
figure(1)
subplot(4,1,1);
plot(xr,mach,\'m\');
title(\'mach number variation along the length of nozzle using non conservative form\');
legend(\'Mach Number\');
grid minor
subplot(4,1,2);
plot(xr,p,\'b\');
title(\'Pressure ratio variation along the length of nozzle using non conservative form\');
legend(\'Pressure ratio\');
grid minor
subplot(4,1,3);
plot(xr,rho,\'g\');
title(\'density ratio variation along the length of nozzle using non conservative form\');
legend(\'Density ratio\');
grid minor
subplot(4,1,4);
plot(xr,t,\'r\');
title(\'Temperature ratio variation along the length of nozzle using non conservative form\');
legend(\'Temperature ratio\');
grid minor
%Conservative Form Solver
g_c=0;
tic
[rho_c,v_c,t_c,g_c,m_c,p_c,mach_c,Mach_throat_c,Pressure_ratio_throat_c,Density_ratio_throat_c,Temperature_ratio_throat_c,k_c]=conservative(n,x,n_throat,dx,a,gamma,g_c,tol,error);
t_consv=toc;
fprintf(\'minimum number of iterations for convergence of conservative form:%f\',g_c);
fprintf(\'time taken by conservstive solver:%f\',t_consv);
%steady state profiles for conservative
figure(2)
subplot(4,1,1);
plot(xr,mach_c,\'m\');
title(\'mach number variation along the length of nozzle using conservative form\');
legend(\'Mach Number\');
grid minor
subplot(4,1,2);
plot(xr,p_c,\'b\');
title(\'Pressure ratio variation along the length of nozzle using conservative form\');
legend(\'Pressure ratio\');
grid minor
subplot(4,1,3);
plot(xr,rho_c,\'g\');
title(\'density ratio variation along the length of nozzle using conservative form\');
legend(\'Density ratio\');
grid minor
subplot(4,1,4);
plot(xr,t_c,\'r\');
title(\'Temperature ratio variation along the length of nozzle using conservative form\');
legend(\'Temperature ratio\');
grid minor
%convergence plots for both the forms
figure(3)
subplot(4,1,1);
plot(k,Mach_throat);
hold on
plot(k_c,Mach_throat_c);
title(\'mach number variation vs iterations using conservative form\');
legend(\'Mach Number Non-Conservative\',\'Mach Number Conservative\');
grid minor
subplot(4,1,2);
plot(k,Pressure_ratio_throat);
hold on
plot(k_c,Pressure_ratio_throat_c);
title(\'Pressure ratio variation vs iteration using conservative form\');
legend(\'Pressure ratio Non Conservative\',\'Pressure ratio Conservative\');
grid minor
subplot(4,1,3);
plot(k,Density_ratio_throat);
hold on
plot(k_c,Density_ratio_throat_c);
title(\'density ratio variation vs iteration using conservative form\');
legend(\'Density ratio Non Conservative\',\'Density ratio Conservative\');
grid minor
subplot(4,1,4);
plot(k,Temperature_ratio_throat);
hold on
plot(k_c,Temperature_ratio_throat_c);
title(\'Temperature ratio variation vs iteration using conservative form\');
legend(\'Temperature ratio Non Conservative\',\'Temperature ratio Conservative\');
grid minor
figure(4)
plot(xr,m);
hold on
plot(xr,m_c);
title(\'Steady Mass Flowrate\');
legend(\'Non-Conservative\',\'Conservative\');
xlabel(\'x/L\');
ylabel(\'Mass Flowrate\');
grid minor
#Non conservative function code
function [rho,v,t,g,m,p,mach,Mach_throat,Pressure_ratio_throat,Density_ratio_throat,Temperature_ratio_throat,k]=non_conservative(n,x,dx,a,gamma,g,tol,error,n_throat)
% calculate initial profiles
rho=1-0.3146*x;
t=1-0.2314*x; % t is temperature
v=(0.1+1.09*x).*t.^0.5;
%time steps
nt=4000;
for k=1:nt
rho_old=rho;
v_old=v;
t_old=t;
m_old=(rho_old.*a.*v_old);
cfl=0.5;
max_vel=v_old+(t_old.^0.5);
dt_vector = (cfl*dx)./(max_vel);
dt = min(abs(dt_vector));
%predictor step
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 eqn
drhodt_p(j)=(-rho(j)*dvdx)-(rho(j)*v(j)*dlogadx)-(v(j)*drhodx);
% momentum eqn
dvdt_p(j)=(-v(j)*dvdx)-(1/gamma)*(dtdx+(t(j)/rho(j))*drhodx);
% energy eqn
dtdt_p(j)=(-v(j)*dtdx)-((gamma-1)*t(j)*(dvdx+(v(j)*dlogadx)));
%soln 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 eqn
drhodt_c(j)=(-rho(j)*dvdx)-(rho(j)*v(j)*dlogadx)-(v(j)*drhodx);
% momentum eqn
dvdt_c(j)=(-v(j)*dvdx)-(1/gamma)*(dtdx+(t(j)/rho(j))*drhodx);
% energy eqn
dtdt_c(j)=(-v(j)*dtdx)-((gamma-1)*t(j)*(dvdx+(v(j)*dlogadx)));
end
% compute avg time derivative
drhodt=0.5*(drhodt_p+drhodt_c);
dvdt=0.5*(dvdt_p+dvdt_c);
dtdt=0.5*(dtdt_p+dtdt_c);
%final soln 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
% apply boundary consitions
%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);
p=rho.*t; %pressure rati
mach=v./(t.^0.5); %mach number
%updated mass flow rate
m=rho.*a.*v;
% number of cycles for convergence
if error>tol
error=max(abs(m-m_old));
g=g+1;
end
% Transient throat values
Mach_throat(g)=mach(n_throat);
Pressure_ratio_throat(g)=p(n_throat);
Density_ratio_throat(g)=rho(n_throat);
Temperature_ratio_throat(g)=t(n_throat);
end
k=linspace(1,g,g);#conservative function code
function [rho,v,t,g_c,m,p,mach,Mach_throat,Pressure_ratio_throat,Density_ratio_throat,Temperature_ratio_throat,k]=conservative(n,x,n_throat,dx,a,gamma,g_c,tol,error)
% Calculate initial profiles
for inip = 1:length(x);
if (x(inip)<=0.5);
rho(inip) = 1.0;
t(inip) = 1.0;
elseif (x(inip)>=0.5) && (x(inip)<=1.5);
rho(inip) = 1.0 - 0.366*(x(inip)-0.5);
t(inip) = 1.0 - 0.167*(x(inip)-0.5);
else
rho(inip) = 0.634 - 0.3879*(x(inip)-1.5);
t(inip) = 0.833 - 0.3507*(x(inip)-1.5);
end
end
v = 0.59./(rho.*a); % velocity
p = rho.*t; % pressure
% Solution vecotors
U1 = rho.*a;
U2 = rho.*a.*v;
U3 = rho.*(t./(gamma-1) + ((gamma/2).*v.^2)).*a;
% Flux vectors
F1 = rho.*a.*v;
F2 = rho.*a.*v.^2 + (1/gamma).*rho.*a;
F3 = rho.*(t./(gamma-1) + (gamma/2).*v.^2).*v.*a + rho.*a.*v;
% Making a copy of all the elements
U1_old = U1;
U2_old = U2;
U3_old = U3;
F1_old = F1;
F2_old = F2;
F3_old = F3;
nt = 4000;
% Outer time loop
for k = 1:nt
% Calculate minium time step
CFL = 0.5;
dt_vector = CFL.*(dx./(t.^0.5 + v));
dt = min(abs(dt_vector));
rho_old = rho;
t_old = t;
v_old = v;
m_old = rho_old.*a.*v_old;
% Making a copy of all the elements
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
dF1dx = (F1(j+1)-F1(j))/dx;
dF2dx = (F2(j+1)-F2(j))/dx;
dF3dx = (F3(j+1)-F3(j))/dx;
dlogadx = (log(a(j+1))-log(a(j)))/dx;
dadx = (a(j+1)-a(j))/dx;
J2 = (1/gamma).*rho(j).*t(j).*dadx;
% Continity Equation
dU1dt_p(j) = -dF1dx;
% Momentum Equation
dU2dt_p(j) = -dF2dx + J2;
% 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
% Primitive variables using solution vector
rho = U1./a; % density
v = U2./U1; % velocity
t = (gamma-1).*((U3./U1) - (gamma/2).*(v.^2)); % temp
p = rho.*t; % pressure
% 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
dF1dx = (F1(j)-F1(j-1))/dx;
dF2dx = (F2(j)-F2(j-1))/dx;
dF3dx = (F3(j)-F3(j-1))/dx;
dlogadx = (log(a(j))-log(a(j-1)))/dx;
dadx = (a(j)-a(j-1))/dx;
J2 = (1/gamma).*rho(j).*t(j).*dadx;
% Continity Equation
dU1dt_c(j) = -dF1dx;
% Momentum Equation
dU2dt_c(j) = -dF2dx + J2;
% Energy Equation
dU3dt_c(j) = -dF3dx;
end
% Compute 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
% Apply boundary conditions
% 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);
% Primitive variables using solution vector
rho = U1./a; % density
v = U2./U1; % velocity
t = (gamma-1).*((U3./U1) - (gamma/2).*(v.^2)); % temp
p = rho.*t; % pressure
% Mass Flow Rate
m = rho.*a.*v;
% Pressure
p = rho.*t; % pressure
% Mach Number
mach = v./(sqrt(t));
if error>tol
error=max(abs(m-m_old));
g_c=g_c+1;
end
% Throat area variables
Mach_throat(g_c)=mach(n_throat);
Pressure_ratio_throat(g_c)=p(n_throat);
Density_ratio_throat(g_c)=rho(n_throat);
Temperature_ratio_throat(g_c)=t(n_throat);
end
k=linspace(1,g_c,g_c);
#solution results
All the parameters listed in table1 and table2 are measured at Throat of nozzle.
Table1: Results of Non conservative form of governing equation and mesh independent study
Table 2: Results of conservative form of governing equation and mesh independent study
Discussion:
1) mesh independent study has been carried out and as you can see from above tables the non conservative form achieves accurate soln for all the parameters listed in table at n=31 itself.
2)from the above tables minimum iterations for convergence increases as the number of grid point increases for both the methods but the increase is more incase of non-conservative form infact at higher grid points the non conservative form require more iterations to convergence compared to conservative form.
3) the non conservative form is much faster compared to conservative form at any number of grid points.
Fig1: steady state plots of non conservative form for different parameters
Fig2: steady state plots of conservative form for different parameters
As you can see from the above figures the steady state plots for both conservative and non conservative form are same for different parameters like mach number, pressure ratio ,density ratio and temperature ratio
Fig3: convergence plots for both conservative and non conservative form for different parameters measured at Throat section of nozzle.
Fig4: steady state plot of non dimensionalized mass flow rate for both conservative and non conservative form of governing equation
as you can see from the above figure the mass flow rate variation over a length of nozzle is much larger for non-conservative form.
case1 n=31; avg mass flow rate for non conservative form is 0.587925.
avg mass flow rate for conservative form is 0.585173.
case2 n=61; avg mass flow rate for non conservative form is 0.583698.
avg mass flow rate for conservative form is 0.58285.
Conclusions
1)conservative form is best in calculating mass flow rate over the length of nozzle compared to non conservative form.
2)non conservative form achieves good accuracy at lower number of grid points only compared to conservative form.so this form is computational viable.
3)non conservative form is faster compared to conservative form, because in conservative form you need to calculate extra variables like flux variables and source variables.
4) as the number of grid point increase the minimum iterations for convergence also increases for both the forms , but the increase in case of non conservative form is more compared to conservative form.( refer above tables)