Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
Aim : Numerical simulation of 1D supersonic nozzle flow using Macormack Method using conservative and non-conservative forms of equation
Conservative Flow :
In this conservative flow , when a fluid is flowing thorugh a pipe a small control volume is taken in some area inside that pipe and this control volume is fixed, so that the equation can be solved to find out the properties of the flow thorugh that point.
Non-Conservative Flow :
In this Non-Conservative flow , a small control volume is taken inside the pipe flow but with a small change, the control volume here is moving with the flow instead of being fixed as in case of conservative flow.
Here we used MacCormack method for solving the equation. This concept is used in programming to give the exact solution.
conventional and nonconventional methods are written in two functions .
Nonconservative Method:
function [v,rho,t,m,p,k,solution_time_1]=nonconservative(n,x,dx,gamma,a,nt,dt)
tic
%area
a=1+2.2*(x-1.5).^2;
nt = 1000;
%calculate initial values
rho=1-0.3146*x;
t=1-0.2314*x;
v=0.1+1.09*x.*t.^0.5;
%outer time loop
for k=1:nt
rho_old = rho;
t_old = t;
v_old = v;
a_old = a;
%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;
%continuitu eq
drhodt_p(j)=-(rho(j)*dvdx)-(rho(j)*v(j)*dlogadx)-(v(j)*drhodx);
%momentum eq
dvdt_p(j)=-(v(j)*dvdx)-((1/gamma)*(dtdx+(((t(j))/rho(j))*drhodx)));
%energy eq
dtdt_p(j)=-(v(j)*dtdx)-((gamma-1)*t(j)*(dvdx+(v(j)*dlogadx)));
%solution update
v(j)=v(j)+(dvdt_p(j)*dt);
rho(j)=rho(j)+(drhodt_p(j)*dt);
t(j)=t(j)+(dtdt_p(j)*dt);
end
%corrector step
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;
%continuitu eq
drhodt_c(j)=-(rho(j)*dvdx)-(rho(j)*v(j)*dlogadx)-(v(j)*drhodx);
%momentum eq
dvdt_c(j)=-(v(j)*dvdx)-((1/gamma)*(dtdx))+((t(j)/rho(j))*drhodx);
%energy eq
dtdt_c(j)=-(v(j)*dtdx)-((gamma-1)*t(j)*(dvdx+(v(j)*dlogadx)));
end
%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 solution update
for i=2:n-1
v(i)=v_old(i)+(dvdt(i)*dt);
rho(i)=rho_old(i)+(drhodt(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));
%mASS FLOW RAtE
m=rho.*a.*v;
%pressure
p = rho.*t;
%for checking convergence of solution
figure(6)
hold on;
if k == 100
plot(x,m,\'g\',\'linewidth\',1.5);
else if k == 200
plot(x,m,\'r\',\'linewidth\',1.5);
else if k == 500
plot(x,m,\'c\',\'linewidth\',1.5);
else if k == 700
plot(x,m,\'b\',\'linewidth\',1.5);
else if k == 1000
plot(x,m,\'y\',\'linewidth\',1.5);
end
end
end
end
end
xlabel(\'x\');
ylabel(\'mass flow rate at different time steps\');
title(sprintf(\'comparision of mass flow rate at different time steps\'));
solution_time_1=toc;
end
end
Conservative Method:
function [v,rho,t,m,p,k,solution_time]=conservative(n,x,dx,gamma,a,nt,dt)
tic
%calculate initial values
for i = 1:n
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);
elseif 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
% Initial condition
v = 0.59./(rho.*a);
U1 = rho.*a;
U2 = (rho.*a).*v;
U3 = rho.*((t./(gamma-1))+ (gamma/2)*v.^2).*a;
v = U2./(rho.*a);
p=rho.*t;
%outer time loop
for k =1:nt
U1_old = U1;
U2_old = U2;
U3_old = U3;
%defining flux value
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 step
for i=2:n-1
J2(i) = (1/gamma).*rho(i).*t(i).*((a(i+1)-a(i))/dx);
%contunity equation
dU1_dt_p(i) = -(F1(i+1)-F1(i))/dx;
%momentum equation
dU2_dt_p(i) = -((F2(i+1)-F2(i))/dx) + J2(i);
%energy equation
dU3_dt_p(i) = -(F3(i+1)-F3(i))/dx;
%solution update
% solution update
U1(i) = U1(i) + dU1_dt_p(i)*dt;
U2(i) = U2(i) + dU2_dt_p(i)*dt;
U3(i) = U3(i) + dU3_dt_p(i)*dt;
end
v=U2./U1;
rho = U1./a;
t = (gamma-1)*((U3./U1) - ((gamma/2)*((U2./U1).^2)));
p=rho.*t;
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 step
for j=2:n-1
J2(j) = (1/gamma)*rho(j).*t(j)*((a(j)-a(j-1))/dx);
dU1_dt_c(j) = -(F1(j)-F1(j-1))/dx;
dU2_dt_c(j) = -((F2(j)-F2(j-1))/dx) + J2(j);
dU3_dt_c(j) = -(F3(j)-F3(j-1))/dx;
end
%avg time derivative
dU1_dt = 0.5*(dU1_dt_p + dU1_dt_c);
dU2_dt = 0.5*(dU2_dt_p + dU2_dt_c);
dU3_dt = 0.5*(dU3_dt_p + dU3_dt_c);
%Final solution update
for i=2:n-1
U1(i) = U1_old(i) + dU1_dt(i)*dt;
U2(i) = U2_old(i) + dU2_dt(i)*dt;
U3(i) = U3_old(i) + dU3_dt(i)*dt;
end
%boundary conditions
%inlet
U1(1) = a(1);
U2(1) = 2*U2(2)-U2(3);
v(1) = U2(1)/U1(1);
t(1) = 1;
U3(1) = U1(1)*((t(1)/(gamma-1)) +((gamma/2)*(v(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);
rho = U1./a;
v = U2./U1;
t = (gamma-1)*((U3./U1) - ((gamma/2)*((U2./U1).^2)));
%MASS FLOW RATE
m=rho.*a.*v;
p=rho.*t;
%for checking convergence of solution
figure(7)
hold on;
if k == 100
plot(x,m,\'g\',\'linewidth\',1.5);
else if k == 200
plot(x,m,\'r\',\'linewidth\',1.5);
else if k == 500
plot(x,m,\'c\',\'linewidth\',1.5);
else if k == 700
plot(x,m,\'b\',\'linewidth\',1.5);
else if k == 1000
plot(x,m,\'y\',\'linewidth\',1.5);
end
end
end
end
end
xlabel(\'x\');
ylabel(\'mass flow rate at different time steps\');
title(sprintf(\'comparision of mass flow rate at different time steps\'));
solution_time=toc;
end
end
Main Program:
clear all
close all
clc
n=31;
x = linspace(0,3,n);
dx = x(2)-x(1);
gamma = 1.4;
%area
a=1+2.2*(x-1.5).^2;
nt = 1000;
dt=0.001;
%non conservative form
[V,R,T,M,P,K,solution_time_1]=nonconservative(n,x,dx,gamma,a,nt,dt);
%conservative form
[v,rho,t,m,p,k,solution_time]=conservative(n,x,dx,gamma,a,nt,dt);
%plots cmparing non conservatve nad consevative parameters
%v plot
figure(1)
plot(x,V,\'r\')
hold on
plot(x,v,\'b\')
xlabel(\'x\');
ylabel(\'velocity\');
title(sprintf(\'comparision of velocity components \'));
legend(\'nonconservative\',\'conservative\');
%rho plot
figure(2)
plot(x,R,\'r\')
hold on
plot(x,rho,\'b\')
xlabel(\'x\');
ylabel(\'density\');
title(sprintf(\'comparision of density components \'));
legend(\'nonconservative\',\'conservative\');
%t plot
figure(3)
plot(x,T,\'r\')
hold on
plot(x,t,\'b\')
xlabel(\'x\');
ylabel(\'temperature\');
title(sprintf(\'comparision of temperature components \'));
legend(\'nonconservative\',\'conservative\');
%m plot
figure(4)
plot(x,M,\'r\')
hold on
plot(x,m,\'b\')
xlabel(\'x\');
ylabel(\'mass flow rate\');
title(sprintf(\'comparision of mass flow rate components \'));
legend(\'nonconservative\',\'conservative\');
%p plot
figure(5)
plot(x,P,\'r\')
hold on
plot(x,p,\'b\')
xlabel(\'x\');
ylabel(\'Pressure\')
title(sprintf(\'comparision of Pressure components \'));
legend(\'nonconservative\',\'conservative\');
Results:
Below are the results for the program.
Plots are drawn by comparing the various parameters for conservative aaaanad nonconservative methods.
Velocity Plot:
Temperature Plot:
Density Plot:
Mass flow rate Plot:
Pressure Plot:
Time Elasped:
Time required for
Conservative time =10.979
Non-Conservative time = 14.945
So, from these above data we can conclude that Non-Conservative form requires more computational time than Conservative form. Conservative form is faster than non conservative form.
The minimum number of cycles a simulation requires to converge depend upon following:
1) Courant Number(C)-which determines the time step size(dt)
dt=min(C*dx/(a+v))
where a=sqrt(T); non dimensionlalized speed of sound
v=velocity; non dimensionalized velocity along the nozzle
2) Error Tolerance/Residual
3) Spacing between two gridpoints(dx)/Number of Gridpoints(n) for a given domain
For the following the minimun number of cycles/minimum number of time steps required by the simulation to run were found out to be the following
From the below plots we can conclude that
Non conservative Form:
From the above graph we can conclude that approximately the solution start converge after the time step 500. As it is shown from the graph the lines of mass of flow start converging after 500 (cyan,blue,yellow lines coincides).
Conservative Form:
From the above graph we can conclude that approximately the solution start converge after the time step 700. As it is shown from the graph the lines of mass of flow start converging after 700 (blue,yellow lines coincides approx).
Grid Dependence Test:
In our case we assumed n=31, when modifying the value of \'n\', the number of time steps should also be modified, so that the elapsed non-dimensional time is same for both cases. Grid independence test can be carried out by comparing the values obtained by conservative and non-conservative forms, for n = 31 & 61. It has to be noted that the non-conservative form results are more accurate when compared with the conservative form results. Even though by doubling the number of grid points, the solutions of the conservative form are slightly closer to the exact value, but still attained the same accuracy of non-conservative form with half as many grid points i.e n = 31.
Flow Parameters - nonconservative - conservative - actual value(from reference table)
n=31 n=61 - n=31 n=61
rho - 0.638 0.636 - 0.645 0.638 - 0.634
t - 0.836 0.835 - 0.839 0.834 - 0.833
p - 0.534 0.530 - 0.538 0.529 - 0.528