Week 7 - Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
Using Macormack Method for simulation of 1-D super-sonic nozzle flow simulation
Conservative Form:- Control Volume is Stationary.
Non-Conservative Form:- Control Volume is dynamic.
One form can be converted into the another form easily and exactly. Conservative form are more accurate but non conserative form are fast.
Governing Equation:-
where V = velocity of the fluid inside the nozzle,
A = Area of the nozzle which varies along the length,
x = length of the nozzle in the x direction ,
rho(p) = density of the fluid and Y = specific heat ratio of the air(= 1.4 assume)
dt = time step
T = Temperature of the air
As we see in the Nozzle that the pressure, velocity ,temperature , area and density varies along the length of the nozzle. In our problem we assume the length of the nozzle of 3 unit and divide into 31 step and with the help of the matlab program we can plot the variation of the variable with the length of the nozzle.
Using the Macormack Method we solve the problem. The macormack method used the predicted step which forward differencing method and the corrected step which is backward differencing method to calculate the value and after averaging the both step the solution is updated after each time step. The forward and backward differencing method is first order approximation method. The central differencing scheme is second order approximation which give more accurate result but it is time consuming. Therefore macormack method is less time consuming method and the result obtained is with less percentage error.
Term used in the conservation form of the matlab program
u = solution vector, F = flux vector, and J = source term
Relation between the solution vector and the flux vector
`(deltaU_1)/(delta t) = -(delta F_1)/(delta x')` - continuity equation
`(deltaU_2)/(delta t) = -(delta F_2)/(delta x') + J_2` -- momentum equation
`(deltaU_3)/(delta t) = -(delta F_3)/(delta x')`-- energy equation
`rho' = U_1/A`-- density of the air
`V' = U_2/U_1`-- velocity of the air
`T' = (gamma -1)*(U_3/U_1 - gamma/2*V^2)`-- temperature of the air
`P = rho'*T'` pressure of the air `
`F_1 = F_1(U_1U_2,U_3)`
`F_2= F_2(U_1U_2,U_3)`
`F_3 = F_3(U_1U_2,U_3)`
Main program
clear all
close all
clc
n = 31; % number of the node point
x = linspace(0,3,n); % the length of the x unit
dx = x(2) - x(1); % the step size of the x unit
% Nonconservation form for n = 31
[rho31_nc, t31_nc, v31_nc, a31_nc, Mass_flow31_nc, Mach_No31_nc, press31_nc, x31_nc,number_iteration_NC]= nonconservation(n)
% conservation form for n = 31
[rho31_c, t31_c, v31_c, a31_c, Mass_flow31_c, Mach_No31_c, press31_c, x31_c,number_iteration_C]= conservation_form(n)
% plot for n = 31
% plot for the non conservation form paramter variable
figure(1)
subplot(3,2,1)
plot(x,rho31_nc,'g')
xlabel('X/L')
ylabel('density non-con')
grid on
subplot(3,2,2)
plot(x,t31_nc,'r')
xlabel('X/L')
ylabel('temperature non-con')
grid on
subplot(3,2,3)
plot(x,v31_nc,'m')
xlabel('X/L')
ylabel('velocity non-con')
grid on
subplot(3,2,4)
plot(x,press31_nc,'b')
xlabel('X/L')
ylabel('pressure non-con')
grid on
subplot(3,2,5)
plot(x,Mass_flow31_nc,'c')
xlabel('X/L')
ylabel('mass rate non-con')
grid on
subplot(3,2,6)
plot(x,Mach_No31_nc,'bl')
xlabel('X/L')
ylabel('mach number non-con')
grid on
% plot the conservation form parameter variable
figure(2)
subplot(3,2,1)
plot(x,rho31_c,'g')
xlabel('X/L')
ylabel('density con')
grid on
subplot(3,2,2)
plot(x,t31_c,'r')
xlabel('X/L')
ylabel('temperature con')
grid on
subplot(3,2,3)
plot(x,v31_c,'m')
xlabel('X/L')
ylabel('velocity con')
grid on
subplot(3,2,4)
plot(x,press31_c,'b')
xlabel('X/L')
ylabel('pressure con')
grid on
subplot(3,2,5)
plot(x,Mass_flow31_c,'c')
xlabel('X/L')
ylabel('mass rate con')
grid on
subplot(3,2,6)
plot(x,Mach_No31_c,'bl')
xlabel('X/L')
ylabel('mach number con')
grid on
figure(3)
plot(x,Mass_flow31_c,'r')
hold on
plot(x,Mass_flow31_nc,'bl')
xlabel('X/L')
ylabel('mass flow variation')
grid on
legend('conservation mass','non-conservation')
title('mass flow variation comparison')
Non conservation form matlab program
function[rho t v a Mass_flow Mach_No press x number_iteration_non_con] = nonconservation(n)
% input data
%n = 31; % number of the node point
x = linspace(0,3,n); % the length of the x unit
dx = x(2) - x(1); % the step size of the x unit
% Initial Condition
rho = 1 - 0.3146*x;
t = 1 - 0.2314*x; % temperature
v = (0.1 +1.09*x).*t.^(0.5);
% area
a = 1+2.2*(x-1.5).^2;
%constant value
gamma = 1.4;
% Time steps
nt = 2000;
c = 0.3;
dt = min((c*dx)./(a+v));
% outer time loop
for i = 1:nt
%dt = min(time_step);
error = 9e9;
tolerance = 1e-4;
number_iteration_non_con = 1;
while(error>tolerance)
% predictor method
rho_old = rho;
t_old = t;
v_old = v;
for j = 2:(n-1)
dvdx = (v(j+1) -v(j))/dx;
drho_dx = (rho(j+1) -rho(j))/dx;
dtdx = (t(j+1)-t(j))/dx;
dlog_a_dx = (log(a(j+1)) - log(a(j)))/dx;
% continuity equation
drhodt_pred(j) = -rho(j)*dvdx -rho(j)*v(j)*dlog_a_dx - v(j)*drho_dx;
%momentum equation
dvdt_pred(j) = -v(j)*dvdx -(1/gamma)*(dtdx + (t(j)/rho(j))*drho_dx);
%energy equation
dtdt_pred(j) = -v(j)*dtdx -(gamma-1)*t(j)*(dvdx +v(j)*dlog_a_dx);
% updated the value
rho(j) = rho(j) + drhodt_pred(j)*dt;
v(j) = v(j) + dvdt_pred(j)*dt;
t(j) = t(j) + dtdt_pred(j)*dt;
end
% for the corrected method
for k = 2 : (n-1)
dvdx = (v(k) -v(k-1))/dx;
drho_dx = (rho(k) -rho(k-1))/dx;
dtdx = (t(k)-t(k-1))/dx;
dlog_a_dx = (log(a(k)) - log(a(k-1)))/dx;
% continuity equation
drhodt_correct(k) = -rho(k)*dvdx -rho(k)*v(k)*dlog_a_dx - v(k)*drho_dx;
%momentum equation
dvdt_correct(k) = -v(k)*dvdx -(1/gamma)*(dtdx + (t(k)/rho(k))*drho_dx);
%energy equation
dtdt_correct(k) = -v(k)*dtdx -(gamma-1)*t(k)*(dvdx +v(k)*dlog_a_dx);
end
% average update value
drhodt_avg = 0.5*(drhodt_pred + drhodt_correct);
dvdt_avg = 0.5*(dvdt_pred + dvdt_correct);
dtdt_avg = 0.5*(dtdt_pred + dtdt_correct);
% for the next time step
for z = 2 : n-1
rho(z) = rho_old(z) + drhodt_avg(z)*dt;
v(z) = v_old(z) + dvdt_avg(z)*dt;
t(z) = t_old(z) + dtdt_avg(z)*dt;
end
error = max(abs(v_old - v));
v_old = v;
number_iteration_non_con = number_iteration_non_con +1;
% apply the boundary condition
% inlet condition
v(1) = 2*v(2) - v(3);
% outlet condition
v(n) = 2*v(n-1) - v(n-2);
t(n) = 2*t(n-1) - t(n-2);
rho(n) = 2*rho(n-1) - rho(n-2);
Mach_No = v./sqrt(t);
Mass_flow = rho.*v.*a;
press = 287*rho.*t;
end
endconservation form matlab program
function [rho, t, v, A, Mass_flow, Mach_No, press, x, number_iteration_con]= conservation_form(n)
% Input data
gamma = 1.4;
n = 31
x = linspace(0,3,n);
dx = x(2) - x(1);
% area of the convergernt-divergent nozzle
A = 1 +2.2.*(x-1.5).^2;
% input parameter
c = 0.5;
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 - 0.366*(x(i)-0.5);
t(i) = 1.0 - 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
% velocity variation in the nozzle throught the domain
v = (0.59)./(rho.*A);
% source vector
U1 = rho.*A;
U2 = rho.*A.*v;
U3 = rho.*((t./(gamma-1)) + (gamma/2)*(v.^2)).*A;
for j = 2:n-1
time_step(j-1) = (c*dx)./(sqrt(t(j)+v(j)));
end
%dt = min(time_step);
dt = 0.0267;
error = 9e9;
tolerance = 1e-4;
number_iteration_con = 1;
while(error>tolerance)
U1_old = U1;
U2_old = U2;
U3_old = U3;
% input flux vector
F1_p = U2;
F2_p = ((U2.^2./U1)+(1-1/gamma)*(U3 - (0.5*gamma*U2.^2./U1)));
F3_p = (gamma*U2.*U3./U1) - ((gamma)*(gamma-1)*0.5*(U2.^3)./(U1.^2));
% for predicted step
for j = 2:(n-1)
dF1dx_p = (F1_p(j+1) - F1_p(j))/dx;
dF2dx_p = (F2_p(j+1)-F2_p(j))/dx;
dJ2dx_p = (1/gamma)*rho(j)*t(j)*(A(j+1) - A(j))/dx;
dF3dx_p = (F3_p(j+1) - F3_p(j))/dx;
% continuity equation
dU1dt_p(j) = -dF1dx_p;
% momentum equation
dU2dt_p(j) = -dF2dx_p + dJ2dx_p;
%energy eqution
dU3dt_p(j) = -dF3dx_p;
% updated the value
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
% corrected step
F1_c = U2;
F2_c = ((U2.^2./U1)+(1-1/gamma)*(U3 - (0.5*gamma*U2.^2./U1)));
F3_c = gamma*(U2.*U3./U1) - (gamma)*(gamma-1)*0.5*((U2.^3)./(U1.^2));
for j = 2:(n-1)
dF1dx_c = (F1_c(j) - F1_c(j-1))/dx;
dF2dx_c = (F2_c(j)-F2_c(j-1))/dx;
dJ2dx_c = (1/gamma)*rho(j)*t(j)*(A(j) - A(j-1))./dx;
dF3dx_c = (F3_c(j) - F3_c(j-1))/dx;
% continuity equation
dU1dt_c(j) = -dF1dx_c;
% momentum equation
dU2dt_c(j) = -dF2dx_c + dJ2dx_c;
%energy eqution
dU3dt_c(j) = -dF3dx_c;
end
% average update value
dU1dt_avg = 0.5*(dU1dt_p +dU1dt_c);
dU2dt_avg = 0.5*(dU2dt_p +dU2dt_c);
dU3dt_avg = 0.5*(dU3dt_p +dU3dt_c);
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
error = max(abs(U1 -U1_old));
U1_old = U1;
number_iteration_con = number_iteration_con+1;
% apply the boundary condition
% inlet 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)/(gamma -1)) + (gamma/2)*(v(1)^2));
% outlet condition
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)- 0.5*gamma*v.^2);
press = rho.*t;
Mass_flow = rho.*A.*v;
Mach_No = v./t.^(0.5);
end
end
Plot for the mass flow rate comparison between the conservation form and non-conservation form.
The mass flow rate variation is more in the non conservative form as compared to the conservative form where mass flow rate more or less constant throught out the length of the nozzle.
Figure 1 The variable inside the nozzle is plotted using the non conservation equation
Figure 2 The variable are plotted using the conservation form of the equation .
As we compare the both figure 1 & figure 2 we can conclude that the pressure, density , velocity and temperature plot are similiar to each other for both the form of the equation i.e conservation form & non conservation form throught the length of the nozzle. But there is a variation in the mass flow rate and the and mach number.
| Using conservation form of the governing equation to obtain the value | |||||
| P/P0 (pressure ratio) | P/P0 (density ratio) | T/T0 (temperature ratio) | Mass flow rate | V/V0 (velocity ratio) | Length Nozzle |
| 1.0291 | 1.0000 | 1.0291 | 0.5958 | 0.1001 | 0 |
| 1.0282 | 0.9999 | 1.0284 | 0.5957 | 0.1121 | 0.1 |
| 1.0274 | 0.9999 | 1.0274 | 0.5955 | 0.1262 | 0.2 |
| 1.0251 | 0.9990 | 1.0262 | 0.5962 | 0.1432 | 0.3 |
| 1.0207 | 0.9964 | 1.0244 | 0.5961 | 0.1634 | 0.4 |
| 1.0137 | 0.9918 | 1.0220 | 0.5959 | 0.1878 | 0.5 |
| 1.0037 | 0.9853 | 1.0187 | 0.5955 | 0.2172 | 0.6 |
| 0.9900 | 0.9764 | 1.0139 | 0.5954 | 0.2532 | 0.7 |
| 0.9713 | 0.9644 | 1.0072 | 0.5952 | 0.2970 | 0.8 |
| 0.9453 | 0.9477 | 0.9975 | 0.5952 | 0.3505 | 0.9 |
| 0.9093 | 0.9245 | 0.9836 | 0.5953 | 0.4155 | 1 |
| 0.8600 | 0.8924 | 0.9637 | 0.5955 | 0.4936 | 1.1 |
| 0.7946 | 0.8491 | 0.9358 | 0.5958 | 0.5857 | 1.2 |
| 0.7118 | 0.7925 | 0.8981 | 0.5960 | 0.6912 | 1.3 |
| 0.6134 | 0.7223 | 0.8492 | 0.5963 | 0.8077 | 1.4 |
| 0.5206 | 0.6523 | 0.7980 | 0.5962 | 0.9140 | 1.5 |
| 0.4040 | 0.5572 | 0.7250 | 0.5962 | 1.0469 | 1.6 |
| 0.3105 | 0.4732 | 0.6562 | 0.5963 | 1.1583 | 1.7 |
| 0.2330 | 0.3957 | 0.5888 | 0.5964 | 1.2580 | 1.8 |
| 0.1726 | 0.3282 | 0.5257 | 0.5966 | 1.3445 | 1.9 |
| 0.1273 | 0.2716 | 0.4689 | 0.5969 | 1.4180 | 2 |
| 0.0943 | 0.2252 | 0.4189 | 0.5971 | 1.4796 | 2.1 |
| 0.0705 | 0.1878 | 0.3755 | 0.5974 | 1.5311 | 2.2 |
| 0.0533 | 0.1577 | 0.3381 | 0.5976 | 1.5741 | 2.3 |
| 0.0408 | 0.1335 | 0.3059 | 0.5979 | 1.6102 | 2.4 |
| 0.0317 | 0.1139 | 0.2783 | 0.5982 | 1.6405 | 2.5 |
| 0.0249 | 0.0981 | 0.2543 | 0.5984 | 1.6665 | 2.6 |
| 0.0199 | 0.0851 | 0.2338 | 0.5988 | 1.6883 | 2.7 |
| 0.0160 | 0.0743 | 0.2154 | 0.5989 | 1.7077 | 2.8 |
| 0.0131 | 0.0655 | 0.2004 | 0.5993 | 1.7233 | 2.9 |
| 0.0107 | 0.0580 | 0.1850 | 0.5997 | 1.7391 | 3 |
| Using Non conservation form of the governing equation to obtain the value | |||||
| P/P0 (pressure ratio) | P/P0 (density ratio) | T/T0 (temperature ratio) | Mass flow rate | V/V0 (velocity ratio) | Length Nozzle |
| 287.0000 | 1.0000 | 1.00000 | 0.0989 | 0.0989 | 0 |
| 286.1296 | 0.9978 | 0.99914 | 0.1120 | 0.1120 | 0.1 |
| 285.7717 | 0.9969 | 0.99877 | 0.1251 | 0.1250 | 0.2 |
| 284.7193 | 0.9943 | 0.99773 | 0.1430 | 0.1428 | 0.3 |
| 283.7688 | 0.9919 | 0.99677 | 0.1626 | 0.1623 | 0.4 |
| 281.9343 | 0.9874 | 0.99493 | 0.1874 | 0.1870 | 0.5 |
| 279.7144 | 0.9818 | 0.99268 | 0.2161 | 0.2153 | 0.6 |
| 276.3491 | 0.9733 | 0.98927 | 0.2521 | 0.2508 | 0.7 |
| 271.9940 | 0.9624 | 0.98478 | 0.2956 | 0.2933 | 0.8 |
| 265.7021 | 0.9464 | 0.97824 | 0.3498 | 0.3460 | 0.9 |
| 257.0005 | 0.9241 | 0.96898 | 0.4159 | 0.4094 | 1 |
| 244.7247 | 0.8923 | 0.95557 | 0.4968 | 0.4857 | 1.1 |
| 228.2284 | 0.8489 | 0.93675 | 0.5944 | 0.5753 | 1.2 |
| 206.9323 | 0.7914 | 0.91102 | 0.7108 | 0.6784 | 1.3 |
| 181.3611 | 0.7202 | 0.87746 | 0.8463 | 0.7927 | 1.4 |
| 153.0133 | 0.6377 | 0.83607 | 0.9997 | 0.9141 | 1.5 |
| 124.4232 | 0.5500 | 0.78820 | 1.1676 | 1.0366 | 1.6 |
| 98.0772 | 0.4640 | 0.73649 | 1.3452 | 1.1544 | 1.7 |
| 75.6420 | 0.3855 | 0.68369 | 1.5279 | 1.2634 | 1.8 |
| 57.6858 | 0.3177 | 0.63267 | 1.7103 | 1.3604 | 1.9 |
| 43.8460 | 0.2613 | 0.58472 | 1.8907 | 1.4458 | 2 |
| 33.4825 | 0.2156 | 0.54124 | 2.0644 | 1.5188 | 2.1 |
| 25.7399 | 0.1787 | 0.50179 | 2.2338 | 1.5824 | 2.2 |
| 20.0380 | 0.1495 | 0.46706 | 2.3937 | 1.6359 | 2.3 |
| 15.7353 | 0.1258 | 0.43568 | 2.5500 | 1.6832 | 2.4 |
| 12.5621 | 0.1072 | 0.40847 | 2.6952 | 1.7225 | 2.5 |
| 10.0828 | 0.0916 | 0.38349 | 2.8397 | 1.7585 | 2.6 |
| 8.2643 | 0.0795 | 0.36218 | 2.9711 | 1.7880 | 2.7 |
| 6.7535 | 0.0688 | 0.34198 | 3.1058 | 1.8162 | 2.8 |
| 5.6881 | 0.0609 | 0.32523 | 3.2245 | 1.8389 | 2.9 |
| 4.6983 | 0.0531 | 0.30847 | 3.3518 | 1.8616 | 3 |
The number of iteration required to converge the solution in the conservation form of the equation is 323 while in the non conservation form of the equation the number of iteration is 2 only. Therefore the non conservation form is fast and quick compared to the conservation form.