Week 7 - Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
Objective
- In this challenge we will simulate the supersonic flow in converging diverging nozzle using Macormack method
- We will use conservative and non conservative governing equations of continuity,momentum and energy
- We will perform Grid Dependence study for conservative and non conservative methods.
Equations for Non Conservative Method
Continuity Equations
`(drho)/dt=-rho*(du)/dx-(rho*u)*(dlogA)/dx- ((u)*(drho)/dx)`
Momentum Equations
`(du)/dt=-(u)*(du)/dx-(1/gamma)*(((dT)/dx)+(T/rho)*(drho)/dx)`
Energy Equation
`(dT)/dt=-(u*(du)/dx)-(gamma-1)*T((du)/dx+(u)*(dlogA)/dx)`
Time Step Calculation
CFL=0.5;
`dt=dx/((sqrt(T)+U)`
Steps for Non Conservative Solution
- Calculate initial variable density area velocity and temperature
- Start with outer loop of iterations 1400
- Differentiate the variables of density velocity and temperature and solve continuity,momentum and energy equations using forward scheme.
- Determine predictor `(du)/dt (dT)/dt (drho)/dt`
- Update density temperature and velocity using `rho=rho_(old)+(du)/dt*(dt)`
- Now use backward scheme and calculate corrector step.`(du)/dt (dT)/dt (drho)/dt`
- Now take average of predictor and corrector and update density temperature and velocity using
- `rho=(rho)_(old)+(((drho)/dt)_(pre)+((drho)/dt)_(co))*0.5*dt`
- Calculate end point velocity using `U(n)=2*U(n-1)-U(n-2)` similarty for Temperature and Densit`
- Find throat at A==1
- Plot mass flow rate as `rho*Area*velocity`
- Calculate `M=(u)/(sqrt T)`
% Supersonic Flow Using Non Conservative Macmorckic Method
clc
clear all
close all
% Prociding Geometrical Cell points
n=31;
x=linspace(0,3,n);
dx=x(2)-x(1);
%Initial the Conditions
rho_initial=1-0.3146*x;
T_initial=1-0.2314*x;
u_initial=(0.1+1.09*x).*T_initial.^(0.5);
% Calculating Areas
A=1+2.2*(x-1.5).^2;
n_throat=find(A-1<0.001);
iter=0;
gamma=1.4;
%Equating Solutions
rho=rho_initial;
T=T_initial;
u=u_initial;
% Time Step Size
cfl=0.5;
for i=1:n
dt_initial(i)=cfl*(dx/(((T_initial(i))^0.5)+(u_initial(i))));
end
nt=1400;
dt = min(dt_initial);
time_required_for_simulation=(nt-1)*dt;
t=0:dt:time_required_for_simulation;
% Discretizing The Govering Equations
for h=1:length(t)
rho_old=rho;
T_old=T;
u_old=u;
for j=2:n-1
% Forward Difference for Predictor step
dudx=(u(j+1)-u(j))/dx;
dAdx=(log(A(j+1))-log(A(j)))/dx;
drhodx=(rho(j+1)-rho(j))/dx;
dTdx=(T(j+1)-T(j))/dx;
% Continuity Equation
drhodt_p(j)=-(rho(j)*dudx)-((rho(j)*u(j))*dAdx)-(u(j)*drhodx);
% Momentum Equation
dudt_p(j)=-(u(j)*dudx)-((1/gamma)*(dTdx+((T(j)/rho(j))*drhodx)));
% Energy Equation
dTdt_p(j)=-(u(j)*dTdx)-((gamma-1)*T(j)*(dudx+(u(j)*dAdx)));
%Modified values of density velocity and temperature
rho(j)=rho_old(j)+(drhodt_p(j)*dt);
T(j)=T_old(j)+(dTdt_p(j)*dt);
u(j)=u_old(j)+(dudt_p(j)*dt);
end
for i=2:n-1
% Backward Difference for Predictor step
dudx=(u(i)-u(i-1))/dx;
dAdx=(log(A(i))-log(A(i-1)))/dx;
drhodx=(rho(i)-rho(i-1))/dx;
dTdx=(T(i)-T(i-1))/dx;
% Continuity Equation
drhodt_c(i)=-(rho(i)*dudx)-((rho(i)*u(i))*dAdx)-(u(i)*drhodx);
% Momentum Equation
dudt_c(i)=-(u(i)*dudx)-((1/gamma)*(dTdx+((T(i)/rho(i))*drhodx)));
% Energy Equation
dTdt_c(i)=-(u(i)*dTdx)-((gamma-1)*T(i)*(dudx+(u(i)*dAdx)));
end
for m=2:n-1
% Calculating Average predictor and corrector eqtations
drhodt_avg(m)=(0.5*dt*(drhodt_p(m)+drhodt_c(m)));
dTdt_avg(m)=(0.5*dt*(dTdt_p(m)+dTdt_c(m)));
dudt_avg(m)=0.5*dt*(dudt_p(m)+dudt_c(m));
%Upadating density velocity and temperatures
rho(m)=rho_old(m)+drhodt_avg(m);
T(m)=T_old(m)+dTdt_avg(m);
u(m)=u_old(m)+dudt_avg(m);
end
u(1) = (2*u(2))-u(3);
% Outlet condition
u(n)=(2*u(n-1))-u(n-2);
rho(n)=(2*rho(n-1))-rho(n-2);
T(n)=(2*T(n-1))-T(n-2);
%calculating
%Mach number
Mach=u./sqrt(T);
%mass flow rate
massflow=rho.*u.*A;
p=rho.*T;
%rho,T,P,M Throat Values
rho_throat(h)=rho(n_throat);
T_throat(h)=T(n_throat);
p_throat(h)=p(n_throat);
M_throat(h)=Mach(n_throat);
% Throat Values
if h==41
plot(x,massflow,'b','linewidth',1.5)
hold on
else if h==111
plot(x,massflow,'r','linewidth',1.5)
hold on
else if h==178
plot(x,massflow,'c','linewidth',1.5)
hold on
else if h==230
plot(x,massflow,'g','linewidth',1.5)
hold on
else if h==450
plot(x,massflow,'*k','linewidth',1.5)
hold on
else if h==750
plot(x,massflow,'.y','linewidth',1.5)
else if h==1400
plot(x,massflow,'v','linewidth',1.5)
hold on
grid on
ylabel('Mass flow Rate','fontweight','bold','fontsize',12);
xlabel('Non dimensional Distance((x))','fontsize',12,'fontweight','bold')
title(sprintf('No of Iterations are=%d',iter))
legend('50 time step','100 time step','150 time step','250 time step','550 time step','700 time step','1400 time step');
end
end
end
end
end
end
end
iter=iter+1;
end
%Exact Solution of Supersonic Flow
rho_exactsol=0.634*ones(1,1400);
p_exactsol=0.528*ones(1,1400);
Mach_exactsol=1*ones(1,1400);
T_exactsol=0.833*ones(1,1400);
%Density
figure(2)
plot(rho_throat,'linewidth',1.5)
hold on
plot(rho_exactsol,'r','linewidth',1.5)
legend('Density ratio simulation','location','northeast','Density ratio Exact','location','northeast');
ylabel('Density Ratio');
xlabel('iterations');
%Temperature
figure(3)
plot(T_throat,'linewidth',1.5)
hold on
plot(T_exactsol,'r','linewidth',1.5)
legend('Temperature ratio simulation','location','northeast','Temperature ratio Exact','location','northeast');
ylabel('Temperature Ratio');
xlabel('iterations');
%Mach Number
figure(4)
plot(M_throat,'linewidth',1.5)
hold on
plot(Mach_exactsol,'r','linewidth',1.5)
legend('Mach Number simulation','location','northeast','Temperature ratio Exact','location','northeast');
ylabel('Mach Number');
xlabel('iterations');
%Pressure
figure(5)
plot(p_throat,'linewidth',1.5)
hold on
plot(p_exactsol,'r','linewidth',1.5)
legend('Pressure Ratio simulation','location','northeast','Pressure ratio Exact','location','northeast');
ylabel('Pressure Ratio');
xlabel('Iterations');
figure(6)
subplot(4,1,1)
plot(x,Mach)
xlabel('Length of nozzle')
ylabel('Mach number','fontsize',10)
grid minor
subplot(4,1,2)
plot(x,p)
xlabel('Length of nozzle')
ylabel('Pressure ratio','fontsize',10)
grid minor
subplot(4,1,3)
plot(x,rho)
xlabel('Length of nozzle')
ylabel('Density ratio','fontsize',10)
grid minor
subplot(4,1,4)
plot(x,T)
xlabel('Length of nozzle')
ylabel('Temp ratio','fontsize',10)
grid minorResults For Non Conservative Methods
- Flow velocity increases as flow transition from Subsonic to Supersonic.
- Pressure reduces and velocity increases in diverging section.
- Density ratio comparison for exact solution and using numerical simulation.
- It can be observe good aggrement between the results.
- Temperature ratio comparison for exact solution and using numerical simulation.
- It can be observe good aggrement between the results.
- Mach Number comparison for exact solution and using numerical simulation.
- It can be observe good aggrement between the results.
- Pressure ratio comparison for exact solution and using numerical simulation.
- It can be observe good aggrement between the results.
- Mass Flow convergence with respect to iterations can be observed below.
Numerical Solution Method for Non conservative approach.
`(du)/dt+(dF)/dx=J`
`u1=rho*A`
`F1=rho*A*u`
`u2=rho*A*u`
`F2=rho*A*U^2+(1/gamma)*p*A`
`(du_(1))/dt=-(dF_(1))/dx`
`(du_(2))/dt=-(dF_(2))/dx+J_(2)`
`(du_(3))/dt=-(dF_(3))/dx`
F1=u2
F2=`u_(2)^2/(u_(1))+(gamma-1)/gamma*(u_(3)-((gamma-1)/2)*(u_(2)^2/(u_(1))))`
`F3=gamma*u_(3)*u_(2)/u_(1) -((gamma-1)/2)*(u_(2)^2/(u_(3)^2)))`
`J2=((gamma-1)/gamma)*(u_(3)- (u_(2)^2)/(2*u_(1)))*(dlogA)/dx`
`rho=u_1/A`
`u=(u_2)/u_1`
`u3=(((T)/(gamma-1))-( gamma/2)*(u_(2)^2))*u_(1)`
`u_(2)=2*u_(2)-u_(3)`
Steps for Conservative Solution
- Calculate initial variable density area velocity and temperature
- Start with outer loop of iterations 1400
- Calculate Fluxes F1,F2,F3 and J2.
- Differentiate the variables velocity using forward scheme.
- Determine predictor `(du)/dt (dT)/dt (drho)/dt`
- Update `u1=u_(old)+(du)/dt*(dt)`
- Now use backward scheme and calculate corrector step.`(du)/dt (dT)/dt (drho)/dt` before that claculate the Flux values F1,F2,F3 and J2.
- Now take average of predictor and corrector and update density temperature and velocity using
- `u1=(u1)_(old)+(((du)/dt)_(pre)+((du)/dt)_(co))*0.5*dt`also for u2 and u3
- Calculate end point velocity using `U(n)=2*U(n-1)-(n-2)` and calculate temperature and density ratio
- Find throat at A==1
- Plot mass flow rate as `rho*Area*velocity`
- Calculate `M=(u)/(sqrt T)`
% Supersonic Flow Using Conservative Macmorckic Method
clc
clear all
close all
n=31;
x=linspace(0,3,n);
dx=x(2)-x(1);
a=1+(2.2*(x-1.5).^2);
n_throat=find(a-1<0.001);
gamma=1.4;
cfl=0.5;
%Total No of Time steps
nt=1400;
for i =1:n
if x(i) <= 0.5
rho_init(i)= 1;
T_init(i) = 1;
elseif x(i) > 0.5 && x(i)<= 1.5
rho_init(i) = 1 - (0.366*(x(i) - 0.5));
T_init(i) = 1 - (0.167*(x(i) -0.5));
elseif x(i) > 1.5 && x(i) <= 3.5
rho_init(i) = 0.634 - (0.3879*(x(i) - 1.5));
T_init(i) = 0.833 - (0.3507*(x(i) -1.5));
end
end
u_init = 0.59./(rho_init.*a);
rho=rho_init;
T=T_init;
u=u_init;
iter=0;
for i=1:n
dt_init(i)=cfl*(dx/(((T_init(i))^0.5)+(u_init(i))));
end
dt = min(dt_init);
time_required_for_simulation=(nt-1)*dt;
t=0:dt:time_required_for_simulation;
% Initial condition
U1_init = rho_init.*a;
U2_init = rho_init.*a.*u_init;
U3_init = rho_init.*a.*((T_init/(gamma-1))+ ((gamma/2)*(u_init.^2)));
U1=U1_init;
U2=U2_init;
U3=U3_init;
for j=1:length(t)
U1_old=U1;
U2_old=U2;
U3_old=U3;
%Fluxes calculation
F1=U2;
F2=((U2.^2)./U1)+(((gamma-1)/gamma)*(U3-((gamma/2)*((U2.^2)./U1))));
F3=(gamma*((U2.*U3)./U1))-(((gamma-1)*gamma)/2)*((U2.^3)./(U1.^2));
for b=2:n-1
dadx=(a(b+1)-a(b))/dx;
J2=(1/gamma)*rho(b)*T(b)*dadx;
%predictor forward differentiation
dU1dt_p(b)=-((F1(b+1)-F1(b))/dx);
dU2dt_p(b)=-((F2(b+1)-F2(b))/dx)+J2;
dU3dt_p(b)=-((F3(b+1)-F3(b))/dx);
%step updating
U1(b)=((dU1dt_p(b)*dt))+U1_old(b);
U2(b)=((dU2dt_p(b)*dt))+U2_old(b);
U3(b)=((dU3dt_p(b)*dt))+U3_old(b);
end
%Inlet
U2(1)=(2*U2(2))-U2(3);
U3(1)=(U1(1))*((T(1)/(gamma-1))+((gamma/2)*u(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);
%variable values
rho=U1./a;
u=U2./U1;
T=(gamma-1)*((U3./U1)-((gamma/2)*u.^2));
%CORRECTOR STEP
F1=U2;
F2=((U2.^2)./U1)+(((gamma-1)/gamma)*(U3-((gamma/2)*((U2.^2)./U1))));
F3=(gamma*((U2.*U3)./U1))-((((gamma-1)*gamma))/2)*((U2.^3)./(U1.^2));
for k=2:n-1
dadx=((a(k))-a(k-1))/dx;
J2=(1/gamma)*rho(k)*T(k)*dadx;
%backward step calculation for corrector loop
dU1dt_c(k)=-((F1(k)-F1(k-1))/dx);
dU2dt_c(k)=-((F2(k)-F2(k-1))/dx)+J2;
dU3dt_c(k)=-((F3(k)-F3(k-1))/dx);
end
for q=2:n-1
%average values for predictor and corrector step
dU1dt_avg(q)=(dU1dt_c(q)+dU1dt_p(q))*0.5*dt;
dU2dt_avg(q)=(dU2dt_c(q)+dU2dt_p(q))*0.5*dt;
dU3dt_avg(q)=(dU3dt_c(q)+dU3dt_p(q))*0.5*dt;
%Final U values Calculation
U1(q)=U1_old(q)+dU1dt_avg(q);
U2(q)=U2_old(q)+dU2dt_avg(q);
U3(q)=U3_old(q)+dU3dt_avg(q);
end
%Inlet
U2(1)=(2*U2(2))-U2(3);
U3(1)=(U1(1))*((T(1)/(gamma-1))+((gamma/2)*(u(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);
%variables
rho=U1./a;
u=U2./U1;
T=(gamma-1)*((U3./U1)-((gamma/2)*u.^2));
p=rho.*T;
Mach=u./(T.^0.5);
massflow=rho.*u.*a;
%value
rho_throat(j)=rho(n_throat);
T_throat(j)=T(n_throat);
p_throat(j)=p(n_throat);
Mach_throat(j)=Mach(n_throat);
iter=iter+1;
if j==60
plot(x,massflow,'b','linewidth',1.5)
hold on
else if j==90
plot(x,massflow,'r','linewidth',1.5)
hold on
else if j==120
plot(x,massflow,'c','linewidth',1.5)
hold on
else if j==240
plot(x,massflow,'g','linewidth',1.5)
hold on
else if j==540
plot(x,massflow,'*k','linewidth',1.5)
hold on
else if j==800
plot(x,massflow,'.y','linewidth',1.5)
hold on
else if j==1400
plot(x,massflow,'*m','linewidth',1.5)
hold on
grid on
ylabel('Mass flow Rate','fontweight','bold','fontsize',12);
xlabel('Non dimensional Distance((x))','fontsize',12,'fontweight','bold')
title(sprintf('No of Iterations are=%d',iter))
legend('50 time step','100 time step','150 time step','250 time step','550 time step','700 time step','1400 time step');
end
end
end
end
end
end
end
end
%Exact Solution of Supersonic Flow
rho_exactsol=0.634*ones(1,1400);
p_exactsol=0.528*ones(1,1400);
Mach_exactsol=1*ones(1,1400);
T_exactsol=0.833*ones(1,1400);
%Density
figure(2)
plot(rho_throat,'linewidth',1.5)
hold on
plot(rho_exactsol,'r','linewidth',1.5)
legend('Density ratio simulation','location','northeast','Density ratio Exact','location','northeast');
ylabel('Density Ratio');
xlabel('iterations');
%Temperature
figure(3)
plot(T_throat,'linewidth',1.5)
hold on
plot(T_exactsol,'r','linewidth',1.5)
legend('Temperature ratio simulation','location','northeast','Temperature ratio Exact','location','northeast');
ylabel('Temperature Ratio');
xlabel('iterations');
%Mach Number
figure(4)
plot(Mach_throat,'linewidth',1.5)
hold on
plot(Mach_exactsol,'r','linewidth',1.5)
legend('Mach Number simulation','location','northeast','Temperature ratio Exact','location','northeast');
ylabel('Mach Number');
xlabel('iterations');
%Pressure
figure(5)
plot(p_throat,'linewidth',1.5)
hold on
plot(p_exactsol,'r','linewidth',1.5)
legend('Pressure Ratio simulation','location','northeast','Pressure ratio Exact','location','northeast');
ylabel('Pressure Ratio');
xlabel('Iterations');
figure(6)
subplot(4,1,1)
plot(x,Mach)
xlabel('Length of nozzle')
ylabel('Mach number','fontsize',10)
grid minor
subplot(4,1,2)
plot(x,p)
xlabel('Length of nozzle')
ylabel('Pressure ratio','fontsize',10)
grid minor
subplot(4,1,3)
plot(x,rho)
xlabel('Length of nozzle')
ylabel('Density ratio','fontsize',10)
grid minor
subplot(4,1,4)
plot(x,T)
xlabel('Length of nozzle')
ylabel('Temp ratio','fontsize',10)
grid minorResults For Conservative Methods
- Flow velocity increases as flow transition from Subsonic to Supersonic.
- Pressure reduces and velocity increases in diverging section.
- Density ratio comparison for exact solution and using numerical simulation.
- It can be observe good aggrement between the results.
- Temperature ratio comparison for exact solution and using numerical simulation.
- It can be observe good aggrement between the results.
- Mach Number comparison for exact solution and using numerical simulation.
- It can be observe good aggrement between the results.
- Pressure ratio comparison for exact solution and using numerical simulation.
- It can be observe good aggrement between the results.
- Pressure ratio comparison for exact solution and using numerical simulation.
- It can be observe good aggrement between the results.
- Mass flow rate vs iterations graph can be mass flow becomes constant before 1400 iterations.
- As we are performed grid independence test the non consrvative part shows higher accuracy than conservative part when comapared with exact solution.
- As we increase the grid points the solution accuracy is improved.
- Line chart shows the non consrvative part has higher accuracy than conservative part when comapared with exact solution.
- As we increase the grid points the solution accuracy is improved.
- Results validated with exact solution.