Week 7 - Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
Objective:
simulate the isentropic flow through a quasi 1D subsonic-supersonic nozzle by deriving conservation and non-conservation forms of the governing equations by implementing MacCormack's technique using Matlab.
Problem statement;
Determine the steady state primitive variables, and analyze the difference between conservative and non conservative form of geverning equations.
Theory
Even though numerical solution of quasi-1-Dimensional, isentropic nozzle flow unreasonable, we are solving this problem to explore conservative and non-conservative form of governing equations and how to approach with MacCormack's technique.
MacCormack's technique is a variant of Lax-Wendroff appraoch but is much simpler in its application. This method is an explicit finate-difference technique which is second-order-accurate in both time and space, this method proceeds by two steps; a predictor step followed by corrected step.
Non conservative Form:
continuty equation: `(partialrho)/(partialt)=(rho)*(partialV)/(partialx)-(rho)*(V)*(partial(lnA))/(partialx)-v*(partialp)/(partialx)`
Momentum Equation : `(partialV)/(partialt)=V*(partialV)/(partialx)-1/(gama )*((partialT)/(partialx)+(T)/(rho)*(partialrho)/(partialx))`
Energy Equation : `(partialT)/(partialt)=V*(partialT)/(partialx)-(gama -1)*T*((partialV)/(partialx)+(V)*(partial(lnA))/(partialx))`
Conservative Form
continuty equation: `(partial(rho*A))/(partialt)+(partial(rho*A*V))/(partialx)=0`
Momentum Equation :`(partial(rho*A*V))/(partialt)+(partial(rho*A*V+(1/(gama)*p*A)))/(partialx)=(1/(gama))*p*((partialA)/(partialx))`
Energy Equation : `(partial(rho((gama)/(gama-1)+(gama)/2V^2)*A))/(partialt)+(partial(rho*A*V+(rho(T/(gama-1)+(gama)/2V^2)V*A)))/(partialx)=0`
Main Program
clear all
close all
clc
L = 3;
n= 31;
x = linspace(0,3,n);
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;
a = 1+2.2*(x-1.5).^2;
[value num] = min(a);
gama = 1.4;
C = 0.5;
for i = 1:n
d_t(i) = C*((dx)/(T(i).^0.5+v(i)));
end
dt = min(d_t);
m = nonconservaive(L,n,x,rho,T,a,v,dt,C,dx,gama,num);
m_1 = nozzelFlow(L,n,x,dx,gama,C);
figure(8)
plot(x,m);
hold on
plot(x,m_1);
xlabel('non dimensional distance through x');
ylabel('mass flow rate');
legend('nonconservative','conservative');non-conservative Form
conservative
function [out] = nozzelFlow(L,n,x,dx,gama,C)
L = 3;
n = 31;
x = linspace(0,3,n);
dx = x(2)-x(1);
gama = 1.4;
C = 0.5;
% initial conditions
for i = 1:length(x)
if x(i) >= 0&& x(i)<=0.5
rho(i) = 1;
T(i) = 1;
a(i) = 1+2.2*(x(i)-1.5).^2;
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);
a(i) = 1+2.2*(x(i)-1.5).^2;
else x(i) >= 1.5&& x(i)<=3
rho(i) = 0.634 - 0.3879*(x(i)-1.5);
T(i) = 0.833-0.3507*(x(i)-1.5);
a(i) = 1+2.2*(x(i)-1.5).^2;
end
end
T = round(T,5);
rho = round(rho,5);
a = round(a,5);
figure(5)
plot(x,a)
[value num] = min(a);
rho_old = rho;
T_old = T;
v = 0.59./(rho_old.*a);
v = round(v,3);
vold = v;
dt = 0.0267;
U1 = round(rho_old.*a,5);
U1_older = U1;
U1_old = U1;
U2 = round(rho_old.*a.*v,5);
U2_old = U2;
U2_older = U2;
U3 = round(rho.*((T./(gama-1))+(gama/2).*(v.^2)).*a,5);
U3_old = U3;
U3_older = U3;
k = 1;
error = 1e09;
mass_old = U2;
nt = 1500;
for k = 1:1600
U1_old = round(rho.*a,5);
U2_old = round(rho.*a.*v,5);
U3_old = round(rho.*((T./(gama-1))+(gama/2).*(v.^2)).*a,5);
k = k ;
U1 = U1_old;
U2 = U2_old;
U3 = U3_old;
F1 = U2;
F2 = round((((U2.^2)./U1)+((gama-1)/gama)*(U3-(gama/2)*((U2.^2)./U1))),5);
F3 = round((((gama.*U2.*U3)./U1) -(gama.*(gama-1)./2).*((U2.^3)./U1.^2)),5);
for i = 2:n-1
J2(i) = ((1/gama).*(rho(i)).*(T(i)).*((a(i+1)-a(i))/dx));
dU1p_dt(i) = -(((F1(i+1)-F1(i))/dx));
dU2p_dt(i) = -((F2(i+1)-F2(i))/dx)+J2(i);
dU3p_dt(i) = -((F3(i+1)-F3(i))/dx);
U1(i) = U1(i)+dU1p_dt(i)*dt;
U2(i) = U2(i)+(dU2p_dt(i))*dt;
U3(i) = U3(i)+(dU3p_dt(i))*dt;
end
rho = round(U1./a,3);
T = round((gama-1)*((U3./U1)-(gama/2)*(U2./U1).^2),5);
v = round(U2./U1,3);
F1 = round(U2,3);
F2 = round((U2.^2)./U1+((gama-1)/gama)*(U3-(gama/2)*((U2.^2)/U1)),5);
F3 = round((gama*U2.*U3)./U1 -(gama*(gama-1)/2)*((U2.^3)./U1.^2),5);
% corrected step
for i = 2:n-1
J2(i) = round((1/gama)*(rho(i))*(T(i))*((a(i)-a(i-1))/dx),5);
dU1c_dt(i) = -round(((F1(i)-F1(i-1))/dx),5);
dU2c_dt(i) = -((F2(i)-F2(i-1))/dx)+J2(i);
dU3c_dt(i) = -round((F3(i)-F3(i-1))/dx,5);
end
dU1_dt = 0.5*(dU1p_dt+dU1c_dt);
dU2_dt = 0.5*(dU2p_dt+dU2c_dt);
dU3_dt = 0.5*(dU3p_dt+dU3c_dt);
for m = 2:n-1
U1(m) = U1_old(m)+((dU1_dt(m))*dt);
U2(m) = U2_old(m)+(dU2_dt(m))*dt;
U3(m) = U3_old(m)+(dU3_dt(m))*dt;
end
% boundery condition
% inlet
U1(1)= a(1);
U2(1)= 2*U2(2)-U2(3);
U3(1)= U1(1)*((T(1)/(gama-1))+(gama/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 = round(U1./a,5);
T = round((gama-1)*((U3./U1)-(gama/2)*(U2./U1).^2),5);
v = round(U2./U1,3);
mass = round(rho.*a.*v,5);
error = max(abs(abs(mass-mass_old)))
mass_old = mass;
v_throat(k) = v(num);
rho_throat(k) = rho(num);
T_throat(k) = T(num);
p_throat(k) = rho(num)*T(num);
mach_throat(k) = v(num)/(T(num))^0.5;
figure(6)
if k == 50
plot(x,mass,'r','linewidth',1.5);
hold on;
elseif k == 1
rho_1 = rho;
v_1 = v;
T_1 = T;
mass_1 = mass;
elseif k == 100
plot(x,mass,'b','linewidth',1.5);
hold on;
elseif k == 200
plot(x,mass,'g','linewidth',1.5);
hold on;
elseif k == 300
plot(x,mass,'m','linewidth',1.5);
elseif k == 500;
plot(x,mass,'y','linewidth',1.5);
hold on;
elseif k == 700
plot(x,mass,'g','linewidth',1.5);
hold on;
elseif k == 1000
plot(x,mass,'m','linewidth',3);
xlabel('non dimensional distance through x');
ylabel('mass flow rate');
legend('100','500','700','1000');
end
end
T_f = T;
v_f = v;
rho = rho;
mass_f = round(rho.*a.*v_f,3);
out = mass_f;
Mach_dista = v_f./(T(num)).^0.5;
p_f = rho.*T_f;
figure(7)
subplot(4,1,1)
plot(x,rho,'r');
xlabel('non dimensional distance');
ylabel('rho/rho_0');
subplot(4,1,2)
plot(x,T_f,'b');
xlabel('non dimensional distance');
ylabel('T/T_0');
subplot(4,1,3)
plot(x,p_f,'b');
xlabel('non dimensional distance');
ylabel('P/P_0');
subplot(4,1,4)
plot(x,Mach_dista,'b');
xlabel('distance');
ylabel('mach number');
figure(8)
subplot(4,1,1)
plot(1:k,rho_throat,'r');
xlabel('time step');
ylabel('rho/rho_0');
subplot(4,1,2)
plot(1:k,T_throat,'b');
xlabel('time step');
ylabel('T/T_0');
subplot(4,1,3)
plot(1:k,p_throat,'b');
xlabel('time step');
ylabel('P/P_0');
subplot(4,1,4)
plot(1:k,mach_throat,'b');
xlabel('time step');
ylabel('mach number');
grid on;
endResult
For conservative :
non dimensionalissed mass flow rate :
Steady state primitives along the length :
primitives variation with time
For non conservative
mass flow rate for different time step
Residuals
Steady state primitives:
primitives variation with times
Mass flow rate calculated by using two diffrerent forms
with 3 decimel accuracy
| `(rho)/(rho o)` | `(T)/(To)` | `(P)/(Po)` | `M` | |
| Analytical result | 0.634 | 0.833 | 0.528 | 1.000 |
| Non-conservative | 0.639 | 0.836 | 0.534 | 0.999 |
| conservative | 0.645 | 0.839 | 0.540 | 0.941 |
Conclusion:
- The effort required for solving conservative form more as compared to non-conservative form.