Simulation of 1D supersonic nozzle using MacCormack method
Objectives
- To write seperate function for conservative and non conservative forms to solve the 1D supersonic nozzle flow equation using Macormack method.
- Figure out the number of cycle for which the simulation must should run for convergence
- compare the normalized flow rate between the conservative and non conservative forms and make suitable conclusion
- Find which method is faster
- Perform the grid dependence test
Theroy
Nozzle:
A nozzle is often a pipe or tube of varying cross sectional area, and it can be used to direct or modify the flow of a fluid (liquid or gas). Nozzles are frequently used to control the rate of flow, speed, direction, mass, shape, and/or the pressure of the stream that emerges from them.
MacCormacks method:
Maccormack method is an explicit finite difference technique which is second order accurate in both time and space .in CFD the maccormack method is widely used discretized scheme for the numerical solution of hyperbolic partial differential equation.
According to the maccormack method.
`v^(t+trianglet).at(i,j)=v^t.at(i,j)+((delv)/(delt))avg.trianglet`
`((delv)/(delt))avg=0.5(((delv)/(delt))predicted+((delv)/(delt))(corrected))`----1
Predictor step:
In the predictor steo the the spacial terms on the right hand side of the eqautionare replaced by tehh forward differnces.
`((delv)/(delt))=((delrho)/(delx))`
`((delv)/(delt))predicted=((rho^tat(i+1,j)-rho^tat(i,j))/(trianglex)-(rho^t(i+1.j)-rho^tat(i,j))/(trianglex))`
`v^(t+trianglet).at(i,j)=v^t.at(i,j)+((delv)/(delt))predicted.trianglet`
Corrector Step:
In the corrector steo, the spatial terms on the R.H.S of the equation are replaaced by backward differentiate
`((delv)/(delt))_(Cerrected)=((rho^tat(i+1,j)-rho^tat(i,j))/(trianglex)-(rho^t(i+1.j)-rho^tat(i,j))/(trianglex))`
After calculating the predictor and corrector step it is put into equation 1
Continuity equation: `(delp)/(delt)=-(rho(delu))/(delx)-(del(lnA))/(delx)-(udelrho)/(delx)`
Momentum equation:
`(delu)/(delt)=-(rho(delu))/(delx)-1/gamma((del(T))/(delx)+(T)/(rho).(delrho)/(delx)`
Energy equation:
`(delT)/(delt)=-(u(delT))/(delx)-(gamma-1)T((del(u))/(delx)+u.(del(ln A))/(delx)`
Conservation form
The governing equation are expressed in a generic form which consist of the solution vector U, the vector F and the source term J.
`(delu)/(delt)+(delF)/(delx)=J`
`U1=rho.A`
`F1=rho.A.u`
`U2=rho.A.u`
`U2=rho.A.u^2+1/gamma.rho.A`
`J2=1/gamma.(delA)/dx`
`U3=rho(e/(gamma-1)+gamma/2.u^2)`
`U3=rho(e/(gamma-1)+gamma/2.u^2)+1/gamma.rho.A`
MAIN CODE
clear all
close all
clc
%define the inputs
n=31;
x=linspace(0,3,n);
nt=1400;
%defining the node of the throat
%16th code
th_node=((n-1)/2)+1;
%calling the conservative approach function
tic
[rho_c,v_c,T_c,mf_c,mac_c,p_c,rho_th_c,T_th_c,v_th_c,M_th_c,mf_th_c,p_th_c]=con_form(x,n,nt,th_node);
simu_con=toc;
%calling the non conservative conservative approach function
tic
[rho_nc,v_nc,T_nc,mf_nc,mac_nc,p_nc,rho_th_nc,T_th_nc,v_th_nc,M_th_nc,mf_th_nc,p_th_nc]=non_con_form(x,n,nt,th_node);
simu_non_con=toc;
%Analytical solution
rho_an=0.634*ones(900,1400);
T_an=0.833*ones(900,1400);
p_an=0.528*ones(900,1400);
M_an=1*ones(900,1400);
%plotting the comparison between the conservative and non
%conservative analyis
figure (7)
plot(x,mf_nc,\'g\')
hold on
plot(x,mf_c,\'r\')
hold on
xlabel(\'Non dimensional distance (x)\')
ylabel(\'Mass flow rate\')
legend(\'NOn conservation\',\'Conservation\')
title_1=sprintf(\'Non conservation analysis time =%f\',simu_non_con)
title_2=sprintf(\'conservation analysis time=%f\',simu_con)
title({\'Comparison of mass flow rate for conservation and conservation approaches\';title_1;title_2})
%ploting the nondimensional parameter between conservative and non conservative]
figure(5)
subplot(5,1,1)
plot(x,mac_nc,\'c\')
hold on
plot(x,mac_c,\'.\',\'color\',\'bla\',\'markersize\',10)
xlabel(\'Non dimensional distance (x)\')
ylabel(\'Mach number (M)\')
legend(\'Non conservation\',\'Conservation\',\'location\',\'northwest\')
title(\'Variation of Non dimenional propertiez along the length of x\')
grid on
subplot(5,1,2)
plot(x,p_nc,\'c\')
hold on
plot(x,p_c,\'.\',\'color\',\'bla\',\'markersize\',10)
xlabel(\'Non dimensional distance (x)\')
ylabel(\'Pressure (P)\')
legend(\'Non conservation\',\'Conservation\',\'location\',\'northwest\')
grid on
subplot(5,1,3)
plot(x,rho_nc,\'c\')
hold on
plot(x,rho_c,\'.\',\'color\',\'bla\',\'markersize\',10)
xlabel(\'Non dimensional distance (x)\')
ylabel(\'Density (rho)\')
legend(\'Non conservation\',\'Conservation\',\'location\',\'northwest\')
grid on
subplot(5,1,4)
plot(x,T_nc,\'c\')
hold on
plot(x,T_c,\'.\',\'color\',\'bla\',\'markersize\',10)
xlabel(\'Non dimensional distance (x)\')
ylabel(\'Temperature ( T )\')
legend(\'Non conservation\',\'Conservation\',\'location\',\'northwest\')
grid on
%plotting the comparison nondimentional parameter at the throat between conservative and non
%conservative analyis
figure(6)
subplot(6,1,1)
plot(linspace(1,nt,nt),M_th_nc,\'r\')
hold on
plot(linspace(1,nt,nt),M_th_c,\'g\')
hold on
plot(linspace(1,nt,nt),M_an,\'b\')
xlabel(\'Time steps\')
ylabel(\'Mach number (M)\')
legend(\'Non conservation\',\'Conservation\')
title(\'Variation of non dimensional flow properities at the throat in accordance with the time steps\')
subplot(6,1,2)
plot(linspace(1,nt,nt),p_th_nc,\'r\')
hold on
plot(linspace(1,nt,nt),p_th_c,\'g\')
hold on
plot(linspace(1,nt,nt),p_an,\'b\')
xlabel(\'Time steps\')
ylabel(\'Pressure(P)\')
legend(\'Non conservation\',\'Conservation\')
grid on
subplot(6,1,3)
plot(linspace(1,nt,nt),rho_th_nc,\'r\')
hold on
plot(linspace(1,nt,nt),rho_th_c,\'g\')
hold on
plot(linspace(1,nt,nt),rho_an,\'b\')
xlabel(\'Time steps\')
ylabel(\'Density(P)\')
legend(\'Non conservation\',\'Conservation\')
grid on
subplot(6,1,4)
plot(linspace(1,nt,nt),T_th_nc,\'r\')
hold on
plot(linspace(1,nt,nt),T_th_c,\'g\')
hold on
plot(linspace(1,nt,nt),T_an,\'b\')
xlabel(\'Time steps\')
ylabel(\'Temperature(T)\')
legend(\'Non conservation\',\'Conservation\')
grid on
Non conservationn function code:
function[rho,v,T,mf,mac,p,rho_th_nc,T_th_nc,v_th_nc,M_th_nc,mf_th_nc,p_th_nc]=non_con_form(x,n,nt,th_node)
%grid spacing
dx=x(2)-x(1);
%initial conditions
rho=1-0.3146*x;
T=1-0.2314*x;
v=(0.1+1.09*x).*(T.^0.5);
gama=1.4;
a=1+(2.2.*((x-1.5).^(2)));
%time loop
for k=1:nt
%time step
d_t_array=0.5.*(dx./((sqrt(T))+v));
dt = min(d_t_array);
%updating the old values at the beginighg
rho_old=rho;
v_old=v ;
T_old=T;
%predictor step
for i=2:n-1
dv_dx=(v(i+1)-v(i))/dx;
dloga_dx=(log(a(i+1))-log(a(i)))/dx;
dT_dx=(T(i+1)-T(i))/dx;
drho_dx=(rho(i+1)-rho(i))/dx;
%continuity equation
drho_dtp(i)=(-rho(i)*dv_dx)-(rho(i)*v(i)*dloga_dx)-(v(i)*drho_dx);
%moemtum
dv_dtp(i)=-v(i)*(dv_dx)-(1/gama)*(dT_dx+(T(i)/rho(i))*drho_dx);
%energy eqaution
dT_dtp(i)=-v(i)*dT_dx-(gama-1)*T(i)*(dv_dx+v(i)*dloga_dx);
%solution update
rho(i)= rho(i)+drho_dtp(i)*dt;
v(i)= v(i)+dv_dtp(i)*dt;
T(i)= T(i)+dT_dtp(i)*dt;
end
%corrector step
for i=2:n-1
dv_dx=(v(i)-v(i-1))/dx;
dloga_dx=(log(a(i))-log(a(i-1)))/dx;
dT_dx=(T(i)-T(i-1))/dx;
drho_dx=(rho(i)-rho(i-1))/dx;
%continuity eqaution
drho_dtc(i)=-(rho(i)*dv_dx)-(rho(i)*v(i)*dloga_dx)-(v(i)*drho_dx);
%momentum eqaution
dv_dtc(i)=-v(i)*dv_dx-(1/gama)*(dT_dx+(T(i)/rho(i))*drho_dx);
%energy equation
dT_dtc(i)=-v(i)*dT_dx-(gama-1)*T(i)*(dv_dx+v(i)*dloga_dx);
end
%calulation of avg time derivative
drho_dt=0.5*(drho_dtp+drho_dtc);
dv_dt=0.5*(dv_dtp+dv_dtc);
dT_dt=0.5*(dT_dtp+dT_dtc);
%update values
for m=2:n-1
rho(m)= rho_old(m)+drho_dt(m)*dt;
v(m)= v_old(m)+dv_dt(m)*dt;
T(m)= T_old(m)+dT_dt(m)*dt;
end
%inlet boundary cond`n
v(1)=2*v(2)-v(3);
%outlet boundary condtion
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);
%mass flow
mf=rho.*a.*v;
%mach number
mac=v./(T.^0.5);
%pressure
p=rho.*T;
%throat values
rho_th_nc(k) = rho(th_node);
v_th_nc(k) = v(th_node);
T_th_nc(k) = T(th_node);
mf_th_nc(k) = mf(th_node);
M_th_nc(k) = mac(th_node);
p_th_nc(k) = p(th_node);
%Residal
rho_res(k) = drho_dt(th_node);
v_res(k) = dv_dt(th_node);
T_res(k) = dT_dt(th_node);
%plotting the mass flow rate at different time time steps along x
figure(2)
if k==50
plot(x,mf,\'b\')
hold on
elseif k==100
plot(x,mf,\'r\')
hold on
elseif k==200
plot(x,mf,\'c\')
hold on
elseif k==550
plot(x,mf,\'g\')
hold on
elseif k==700
plot(x,mf,\'y\')
hold on
elseif k==1000
plot(x,mf,\'m\')
hold on
elseif k==1400
plot(x,mf,\'--\',\'color\',\'bla\')
hold on
grid on
xlabel(\'Non dimensional distance through the nozzle x \\rightarrow\')
ylabel(\'Non dimensional mass flow rate x \\rightarrow\')
legend (\'50\\deltat\',\'100\\deltat\',\'200\\deltat\',\'550\\deltat\',\'700\\deltat\',\'1000\\deltat\',\'1400\\deltat\')
title(sprintf(\' Mass flow throught the nozzle at different time step Non conservative approach \'))
end
end
%plotting the residual along the time step
figure(3)
plot(linspace(1,nt,nt),rho_res,\'r\')
hold on
plot(linspace(1,nt,nt),v_res,\'g\')
hold on
plot(linspace(1,nt,nt),T_res,\'b\')
grid on
xlabel(\'Time steps\')
ylabel(\'Residual\')
legend(\'Density plot \',\'velocity plot\',\'Temperature plot\')
title(\'Residual plot for Non conservation form\')
endConservation code:
function [rho,v,T,mf,mac,p,rho_th_c,T_th_c,v_th_c,M_th_c,mf_th_c,p_th_c]=con_form(x,n,nt,th_node)
% grid spacing
dx=x(2)-x(1);
gama=1.4;
%intial profile
for i= 1:n
if (x(i)>=0 && x(i)<=0.5)
rho(i) = 1;
T(i) =1;
end
if (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);
end
if (x(i)>=1.5 && x(i)<=3.5)
rho(i) = 0.634-0.3879*(x(i)-1.5);
T(i) =0.833-0.3507*(x(i)-1.5);
end
end
%inital condition
a=1+(2.2.*((x-1.5).^(2)));
v=0.59./(rho.*a);
u1=rho.*a;
u2=rho.*a.*v;
u3=u1.*((T./(gama-1))+((gama/2)*(v.^2)));
%Time loop
for k=1:nt
CFL_no =0.5;
d_t=min(CFL_no.*(dx./((T.^0.5)+v)));
%time loop
u1_old=u1;
u2_old=u2;
u3_old=u3;
%flux term
f1=u2;
term=(((gama-1)/gama)*(u3-((gama/2)*((u2.^2)./u1))));
f2=((u2.^2)./u1)+term;
f3=((gama*u2.*u3)./u1)-(gama*(gama-1)*(u2.^3))./(2*u1.^2);
%predictor step
for j=2:n-1
df1dx(j)=(f1(j+1)-f1(j))/dx;
df2dx(j)=(f2(j+1)-f2(j))/dx;
df3dx(j)=(f3(j+1)-f3(j))/dx;
dadx(j)=(a(j+1)-a(j))/dx;
j2(j)=(1/gama).*rho(j).*T(j).*dadx(j);
%continuity step\'du1dt_
du1_dtp(j)=-df1dx(j);
%momentum eqaution
du2_dtp(j)=-df2dx(j)+j2(j);
%energy eqaution
du3_dtp(j)=-df3dx(j);
%solution update
u1(j)=u1(j)+du1_dtp(j)*d_t;
u2(j)=u2(j)+du2_dtp(j)*d_t;
u3(j)=u3(j)+du3_dtp(j)*d_t;
end
%primitives
rho=u1./a;
v=u2./u1;
T=(gama-1)*((u3./u1)-((gama/2)*(v.^2)));
%flux update
f1=u2;
term=(((gama-1)/gama)*(u3-((gama/2)*((u2.^2)./u1))));
f2=((u2.^2)./u1)+term;
f3=((gama*u2.*u3)./u1)-(gama*(gama-1)*(u2.^3))./(2*u1.^2);
%corrector step
for m=2:n-1
df1dx(m)=(f1(m)-f1(m-1))/dx;
df2dx(m)=(f2(m)-f2(m-1))/dx;
df3dx(m)=(f3(m)-f3(m-1))/dx;
dadx_1(m)=(a(m)-a(m-1))/dx;
j2(m)=(1/gama).*rho(m).*T(m).*dadx_1(m);
%contiuuity eqaution
du1_dtc(m)=-df1dx(m);
%momentum eqaution
du2_dtc(m)=-df2dx(m)+j2(m);
%energy eqaution
du3_dtc(m)=-df3dx(m);
end
%avergae values
du1dt=0.5*(du1_dtp+du1_dtc);
du2dt=0.5*(du2_dtp+du2_dtc);
du3dt=0.5*(du3_dtp+du3_dtc);
%final update
for p =2:n-1
u1(p)=u1_old(p)+du1dt(p)*d_t;
u2(p)=u2_old(p)+du2dt(p)*d_t;
u3(p)=u3_old(p)+du3dt(p)*d_t;
end
%inlet boundary condition
u1(1)=rho(1).*a(1);
u2(1)=2*u2(2)-u2(3);
v(1)=u2(1)/u1(1);
u3(1)=u1(1).*((T(1)./(gama-1))+((gama/2)*(v(1).^2)));
%outflow boundary conditon
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);
%primitives
rho=u1./a;
v=u2./u1;
T=(gama-1)*((u3./u1)-((gama/2)*(v.^2)));
%updates
mf=rho.*v.*a;
mac=v./(sqrt(T));
p=rho.*T;
%Throat values
rho_th_c(k) = rho(th_node);
v_th_c(k) = v(th_node);
T_th_c(k) = T(th_node);
mf_th_c(k) = mf(th_node);
M_th_c(k) = mac(th_node);
p_th_c(k) = p(th_node);
%Residal
u1_res(k) = du1dt(th_node);
u2_res(k) = du2dt(th_node);
u3_res(k) = du3dt(th_node);
%plotting the mass flow rate at different time time steps along x
figure(1)
if k==50
plot(x,mf,\'b\')
hold on
elseif k==100
plot(x,mf,\'r\')
hold on
elseif k==200
plot(x,mf,\'c\')
hold on
elseif k==550
plot(x,mf,\'g\')
hold on
elseif k==700
plot(x,mf,\'y\')
hold on
elseif k==1000
plot(x,mf,\'m\')
hold on
elseif k==1400
plot(x,mf,\'--\',\'color\',\'bla\')
hold on
grid on
xlabel(\'Non dimensional distance through the nozzle x \\rightarrow\')
ylabel(\'Non dimensional mass flow rate x \\rightarrow\')
legend (\'50\\deltat\',\'100\\deltat\',\'200\\deltat\',\'550\\deltat\',\'700\\deltat\',\'1000\\deltat\',\'1400\\deltat\')
title(sprintf(\' Mass flow throught the nozzle at different time step conservative approach \'))
end
end
%plotting the residual
figure(4)
plot(linspace(1,nt,nt),u1_res,\'r\')
hold on
plot(linspace(1,nt,nt),u2_res,\'g\')
hold on
plot(linspace(1,nt,nt),u3_res,\'b\')
grid on
xlabel(\'Time steps\')
ylabel(\'Residual\')
legend(\'u1 plot \',\'u2 plot\',\'u3 plot\')
title(\'Residual plot for conservation form\')
end
PLOTS:Non Conservation plot for mass flow rate at 61 ānd 31 node along the lenght of the nozzle at different time step
At node:31
At node :61
From the above two plots of conservation approach at 61 and 31 node along the lenght of the nozzle we can see the first line that is at 50 time step is way beyond convergence but as we proceede through increase in time steps we conclude that at 700 time step represented by yellow line, the solution approximatly gets converged and secondly it can also be stated it is the minimum number of cycle needed to get the solution get converged.
PLOTS:Conservation plot for mass flow rate at 61ānd 31 node along the lenght of the nozzle at different time step
NOTE:31
AT NODE:61
Here the most of the lines are close to each other and thus aprroaching to the convergence of unlike that of the non conservative approach.The 700 time step represented bt yellow colour give the first evidence that the solution is approaching to convergence and if we proceede the next timestep after 700 time ie 1000 and 1400 steps the lines coincide with each other.
PLOT:Mass flow rate plot between conservation and non-conservation plot
31nodes
61 nodes
The non-conservation plots looks a lot more varying in contrast to consservation wich have nearly the smoth lines.Another thing to be noticed is that at the throat value the mass flow decreased and then increases it as it crosses.The reason behind these two point is the that the control point is also moving with flow therby resulting in the variation in parameter as the flow advances but on the other side in conservative approach the control is fixed thus mainting the near to constant mass flow rate.This contribute to showcasing a smooth line
From the above we can also seet the the analysis time for both conservative and non conservative approaches.we see that non conservatives converges the solution at faster rate with 3.204sec with 61 nodes and 3.45 with 31 nodes while for the comservative appraoch the time taken to converge are 3.6 and4.0.12 with 31 and 61 nodes respectively.The factors contributing to the longer time to converge in conservative approach is additional terms such as the flux term and solution term which are required to decode the primitive terms.
With 31 nodes
with 61 nodes
The above plots are of the variation of paramter along the length the of the nozzle x.The parameter taken are Mach number,Pressure,Density,Temperature.The profile for both conservative and non conservative are same and matches with the condition of convergent and divergent nozzle.Talking the about the mach number is it gives a perfection variation of an actual nozzle with initally being low at the convergent side and then the profile increase as it advances towards towards the divergent sise of nozzle which in the supersonic state.The pressure,density and temperature are high at the inital but keeps on decreasing at the end of the nozzle.
Variation of parameters at the throat of the nozzle
The profile of the four flow properties for Non conservative approach looks a variying initally and soon the it becomes stable and smooth at it reaches 210 time step.in the contrast to that of Non conservative, the Conservative approaches at the starting of the of the time step shos a fairly smooth profile for all the flow properties and eventually both the plots coincides
Grid dependence test
The grid dependence test is neccessary it show how the relation of solution is with number of grids.For this we have computed with 31 grids and 61 grids the came the results shows that the solution are fairly same and for this plots are are taken with both number of grids.This conclude the fact that the relation of solution with number of grid point is independent.
Non conservative
No of grids Density Temperature Pressure Mach number
31 0.638 0.8365 0.534 0.999
61 0.637 0.8354 0.532 0.999
Conservative
31 0.365 0.840 0.546 0.989
61 0.638 0.835 0.533 0.999