Week 7 - Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
Introduction to Nozzle
A nozzle is a device designed to control the direction or characteristics of a fluid flow (to increase velocity) as it exists (or enters) an enclosed chamber or pipe. 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 pressure of the stream that emerges from them. In a nozzle, the velocity of the fluid increases at the expense of its pressure energy.
Flow through nozzles occurs in a variety of engineering applications:
- Convergent nozzles are used for subsonic and sonic flows. They can also be used as flow measuring and flow regulating devices.
- Convergent-Divergent nozzles are used for Supersonic flows. It is comparitively difficult to obtain convergent-divergent passages in compressors and turbine blade rows.
In this project Convergent-Divergent nozzle is concidered code is written to solve the 1D Supersonic nozzle flow equations using Macormack Method.
shows the flow from an infinite reservoir through a convergent-divergent nozzle to an exhaust chamber. Though the main subject here is isentropic flow,it will be observed that flow through a convergent-divergent nozzle is not isentropic under all condition.The static pressure distributions for various values of the pressure ratio pbp0">pbp0pbp0 are shown in curves a,b,c,...,h,i and j. In curves a and b the pressure ratio pbp0">pbp0pbp0,across the nozzle is such that the flow is accelerating only up to the throat;the diverging part acts as a diffuser through which the pressure rises to the back pressure pb">pbpb. In curve critical conditions(pt=p;Mt=1">pt=p;Mt=1pt=p;Mt=1) are reached at the throat but the back pressure is such that the diverging part still acts as a diffuser.
Curve h corresponds to the design value of the back pressure.All other curves are for off-design conditions.Here the flow is supersonic in the entire divergent part of the nozzle.The flow in the curves a,b,c,...,h,i and j is isentropic but non-isentropic in other curves. After the conditions of curve c, when the back pressure is further lowered expansion takes
place to supersonic velocity beyond the throat to a point where a discontinuity in the flow occurs; this is on account of the"off design”"value of the back pressure and the supersonic flow as in curves d and e.In these cases at some point the minimum pressure reached in supersonic expansion has to rise to the back pressure. But it is known that for the
deceleration of supersonic flow the passage should be convergent. Therefore, the existing shape of the nozzle downstream of this point is incompatible with the required process;as a result of this the flow readjusts itself to the shape of the passage by suddenly becoming subsonic. The pressure,temperature and density suddenly rise to values which are compatible with the subsonic flow immediately downstream of the supersonic flow region.Such a sudden change of supersonic flow to subsonic occurs through a plane of discontinuity between the supersonic and subsonic flow regions.This plane of discontinuity is called a shock wave. The flow through the shock wave is no longer isentropic.
When the back pressure is lowered further the shock wave moves downstream till it reaches the exit as in curves f and g;here the back pressure is still higher than the design value.The nozzle exit pressure pe">pepe,rises to the back pressure pb">pbpb,through a shock wave outside the nozzle exit section as shown in curves f and g. Curves i and j correspond to back pressure lower than the nozzle exit pressure.The conditions in the nozzle remain unchanged and the adjustment of the nozzle exit pressure to the exhaust chamber pressure occurs through an expansion wave outside the nozzle.
The Maccormack method-
This method is employed to solve the governing equations of quasi 1D supersonic nozzle flow rather than the differencing schemes as we have studied earlier. Maccormack method is an explicit finite difference technique which is second order accurate in both time and space.In computational fluid dynamics, the MacCormack method is a widely used discretization scheme for the numerical solution of hyperbolic partial differential equations. This method essentially uses a predictor step and a corrector step which will be depicted in the code below.
According to the maccormack method, the value of variable at grid point i,j at time level t+dt is given as follows-
vi,jt+Δt=vi,jt+(∂v∂t)avg⋅Δt">vt+Δti,j=vti,j+(∂v∂t)avg⋅Δtvi,jt+Δt=vi,jt+(∂v∂t)avg⋅Δt
where the average time derivative is the mean value between time t and t+dt. The average time derivative is obtained from the predictor-corrector philosophy.
(∂v∂t)avg=0.5[(∂v∂t)predected+(∂v∂t)corected]">(∂v∂t)avg=0.5[(∂v∂t)predected+(∂v∂t)corected](∂v∂t)avg=0.5[(∂v∂t)predected+(∂v∂t)corected]
Predictor Step-
In the predictor step, the spatial terms on the right-hand side of the equation are replaced by forward differences. The RHS of the equation is known as all the values of variables are for time-step 't'. Next, the value of the variable on the LHS is calculated by employing the time marching solution and for time-step 't+dt' value of the variable at all grid points is obtained.
Consider the following equation for example-
(∂v∂t)=(∂ρ∂x)">(∂v∂t)=(∂ρ∂x)(∂v∂t)=(∂ρ∂x)
(∂v∂t)predected=(ρi+1,jt-ρi,jtΔx+ρi+1,jt-ρi,jtΔy)">(∂v∂t)predected=(ρti+1,j−ρti,jΔx+ρti+1,j−ρti,jΔy)(∂v∂t)predected=(ρi+1,jt-ρi,jtΔx+ρi+1,jt-ρi,jtΔy)
vi,jt+Δt=vi,jt+(∂v∂t)predicted⋅Δt">vt+Δti,j=vti,j+(∂v∂t)predicted⋅Δtvi,jt+Δt=vi,jt+(∂v∂t)predicted⋅Δt
This is the predicted value of variable which will be used in the corrector step to evaluate the corrected value of variable.
Corrector Step-
In the corrector step, the spatial terms on the right-hand side are replaced by backward differences. The value of variables in RHS are taken which are obtained in the predicted step for time-step 't+dt'. Then the corrected time derivative is calculated.
(∂v∂t)corrected=(ρi,jt+Δt-ρi-1,jt+ΔtΔx+ρi,jt+Δt-ρi-1,jt+ΔtΔy)">(∂v∂t)corrected=⎛⎝ρt+Δti,j−ρt+Δti−1,jΔx+ρt+Δti,j−ρt+Δti−1,jΔy⎞⎠(∂v∂t)corrected=(ρi,jt+Δt-ρi-1,jt+ΔtΔx+ρi,jt+Δt-ρi-1,jt+ΔtΔy)
Now, the predicted and corrected value of time derivative is used to calculate the average time derivative which in turn gives the value of the variable at a given point for a particular time-step. The same steps are repeated until the values of variables are obtained for all grid points for a given number of time-steps
In this project the fluid properties are studied as it passes through the convergent-divergent nozzle. The fluid properties are studied for both Conservative and Non-Conservative types of flows. Though both conservative and non-conservative are mathematically equal, they are different in numerical descritisation techniques. One form can be converted to another exactly.
Governing equations for Quasi 1-D Supersonic nozzle flow-
the governing contunity,momentum and energy equations for Quasi 1D steady icentropic flow is given by
Contunity Equation: ρ1V1A1=ρ2V2A2">ρ1V1A1=ρ2V2A2ρ1V1A1=ρ2V2A2
Momentum Equation: ρ1A1+ρ1V12+∫pdV=p2A2+ρ1V22">ρ1A1+ρ1V21+∫pdV=p2A2+ρ1V22ρ1A1+ρ1V12+∫pdV=p2A2+ρ1V22
Energy Equation: h1+v122=h2+v222">h1+v212=h2+v222h1+v122=h2+v222
The conservative and non-conservative form of governing equations are derived using the integral form of continuity, momentum and energy equation.The derivation is given in the book- Computational Fluid Dynamics by John Anderson Jr.
Non- Conservative Form
Contunity Equation: ∂ρ∂t=-ρ∂u∂x-ρu∂(lnA)∂x-u∂ρ∂x">∂ρ∂t=−ρ∂u∂x−ρu∂(lnA)∂x−u∂ρ∂x∂ρ∂t=-ρ∂u∂x-ρu∂(lnA)∂x-u∂ρ∂x
Momentum Equation: ∂u∂t=-u∂u∂x-1γ(∂T∂x+Tρ∂ρ∂x)">∂u∂t=−u∂u∂x−1γ(∂T∂x+Tρ∂ρ∂x)∂u∂t=-u∂u∂x-1γ(∂T∂x+Tρ∂ρ∂x)
Energy Equation:∂T∂t=-u∂T∂x-(γ-1)T(∂u∂x+u∂(lnA)∂x)">∂T∂t=−u∂T∂x−(γ−1)T(∂u∂x+u∂(lnA)∂x)∂T∂t=-u∂T∂x-(γ-1)T(∂u∂x+u∂(lnA)∂x)
All the variables are the non-dimensional values. These equations are solved using Maccormack method and the steady-state values of all variables i.e. Temperature, Velocity and density are obtained. From these variables, the Mach number and pressure is obtained.
p=ρT">p=ρTp=ρT
Mach number = u/T">√TT
Calculation of Time-step-
In order for the solution to be stable, the concept of Courant flow number was introduced in previous projects. For a linear Hyperbolic equation, the value of CFL<1 will give a stable solution.
For this project, CFL=0.5 and value of time-step is calculated as follows-
Δt=CΔxV+a">Δt=CΔxV+aΔt=CΔxV+a
where, C is the Courant number, V is the local flow velocity at a point and a is the local speed of sound. Value of dt at all grid points is calculated and the minimum of them is chosen for the simulation.
Boundary Conditions-
At the inlet:
Velocity at the inlet is changing with time, while the density and temperature are constant with respect to time.
u1=2u2+u3">u1=2u2+u3u1=2u2+u3
Velocity at the inlet is calculated by the method of linear extrapolation. Values of velocity at node 2 and 3 are used for the same.
At the Outlet:
All the variables at the outlet are kept changing with respect to time and linear extrapolation is used to calculate the variables at the outlet.
uN=2uN-1+uN-2">uN=2uN−1+uN−2uN=2uN-1+uN-2
TN=2TN-1+TN-2">TN=2TN−1+TN−2TN=2TN-1+TN-2
ρN=2ρN-1+ρN-2">ρN=2ρN−1+ρN−2ρN=2ρN-1+ρN-2
Initial Conditions-
Area of the nozzle will remain constant throughout the time of simulation.
A=1+2.2(x-1.5)2">A=1+2.2(x−1.5)2A=1+2.2(x-1.5)2
The flow-field variables vary with respect to time. At t=0, their values are given as-
ρ1+0.3146x">ρ1+0.3146xρ1+0.3146x
T=1+0.2314x">T=1+0.2314xT=1+0.2314x
u=(0.1+1.09x)T0.5">u=(0.1+1.09x)T0.5u=(0.1+1.09x)T0.5
Conservative Form-
The governing equations are expressed in a generic form which consists of the solution vector U, the flux vector F and the source term J.
∂U∂t+∂F∂x=J">∂U∂t+∂F∂x=J∂U∂t+∂F∂x=J
The equations for Quasi 1-D supersonic flow can be expressed in a similar generic form as follows-
U1=ρA">U1=ρAU1=ρA F1=ρAu">F1=ρAuF1=ρAu
U2=ρAu">U2=ρAuU2=ρAu F2=ρAu2+1γpA">F2=ρAu2+1γpAF2=ρAu2+1γpA J2=1γp∂A∂x">J2=1γp∂A∂xJ2=1γp∂A∂x
U3=ρ(eγ-1+γ2u2)">U3=ρ(eγ−1+γ2u2)U3=ρ(eγ-1+γ2u2) F3=ρ(eγ-1+γ2u2)+1γpA">F3=ρ(eγ−1+γ2u2)+1γpAF3=ρ(eγ-1+γ2u2)+1γpA
The equations which are solved using Maccormack method are-
∂U1∂t=∂F1∂x">∂U1∂t=∂F1∂x∂U1∂t=∂F1∂x ∂U2∂t=-∂F2∂x+J2">∂U2∂t=−∂F2∂x+J2∂U2∂t=-∂F2∂x+J2 ∂U3∂t=-∂F3∂x">∂U3∂t=−∂F3∂x∂U3∂t=-∂F3∂x
Further, the flux vectors and source terms are expressed in terms of solution vectors. Also, to find the field variables, they are written in terms of solution vectors.
F1=U1">F1=U1F1=U1
F2=U22U1+γ-1γ(U3-γ2U22U1)">F2=U22U1+γ−1γ(U3−γ2U22U1)F2=U22U1+γ-1γ(U3-γ2U22U1) F3=γU2U3U1+γ(γ-1)2U23U12">F3=γU2U3U1+γ(γ−1)2U32U21F3=γU2U3U1+γ(γ-1)2U23U12 J2=γ(γ-1)γ(U3-γ2U22U1)∂lnA∂x">J2=γ(γ−1)γ(U3−γ2U22U1)∂lnA∂xJ2=γ(γ-1)γ(U3-γ2U22U1)∂lnA∂x
Also,
ρ=U1A">ρ=U1Aρ=U1A u=U2U1">u=U2U1u=U2U1 T=(γ-1)(U3U1-γ2u2)">T=(γ−1)(U3U1−γ2u2)T=(γ-1)(U3U1-γ2u2)
Calculation of time step is the same as that of the procedure followed for non-conservation form.
Boundary conditions-
At inlet:
U2 and U3 are changing with respect to time whereas U1 is kept as constant.
U2(i=1)=U2(i=2)-U2(i=3)">U2(i=1)=U2(i=2)−U2(i=3)U2(i=1)=U2(i=2)-U2(i=3)
U3=U1(Tγ-1+γ2u2)">U3=U1(Tγ−1+γ2u2)U3=U1(Tγ-1+γ2u2)
At outlet:
U1(n)=2U1(n-1)-U1(n-2)">U1(n)=2U1(n−1)−U1(n−2)U1(n)=2U1(n-1)-U1(n-2)
U2(n)=2U2(n-1)-U2(n-2)">U2(n)=2U2(n−1)−U2(n−2)U2(n)=2U2(n-1)-U2(n-2)
U3(n)=2U3(n-1)-U3(n-2)">U3(n)=2U3(n−1)−U3(n−2)U3(n)=2U3(n-1)-U3(n-2)
Initial conditions-
Area of the nozzle is constant throughout the runtime of simulation and its value is same as that taken for non-conservative form.
Initial velocity is given as,
u=0.59ρA">u=0.59ρAu=0.59ρA
At t=0, values of density and temperature throughout the nozzle are given as follows,
">rho'=1.0">T'=1.0` for 0≤x′≤0.5">0≤x'≤0.50≤x′≤0.5
ρ′=1.0-0.366(x′-0.5)">ρ'=1.0−0.366(x'−0.5)ρ′=1.0-0.366(x′-0.5) T′=1.0-0.167(x′-0.5)">T'=1.0−0.167(x'−0.5)T′=1.0-0.167(x′-0.5) for 0.5≤x′≤1.5">0.5≤x'≤1.50.5≤x′≤1.5
ρ′=0.634-0.387(x′-1.5)">ρ'=0.634−0.387(x'−1.5)ρ′=0.634-0.387(x′-1.5) T′=0.833-0.3507(x′-1.5)">T'=0.833−0.3507(x'−1.5)T′=0.833-0.3507(x′-1.5) for 1.5≤x′≤3.5">1.5≤x'≤3.5
Matlab Code
For Non-conservative Approach:
function[v,rho,T,mass,M,p,v_throat,rho_throat,T_throat,mass_throat,M_throat,p_throat]=non_conservative_approach(n,L,cfl,nt,th)
%inputs
n=31;
x=linspace(0,3,n);
dx=x(2)-x(1);
%conidtions
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;
th=find(a==1);
gama=1.4;
%time steps
%Time Step calculation based on Courant number.
cfl=0.5;
dt = min(cfl.*(dx./((sqrt(T))+v)));
%Number of iterations
nt=1400;
%rho_initial=rho;
%T_initial=T;
%v_initial=v;
for k=1:nt
%copy of the variables
rho_old=rho;
T_old=T;
v_old=v;
%predictor step
for j= 2:n-1
dv_dx= (v(j+1)-v(j))/dx;
drho_dx =(rho(j+1)-rho(j))/dx;
dTdx = (T(j+1)-T(j))/dx;
da_dx =(log(a(j+1))-log(a(j)))/dx;
%continuity
drhodt_p(j) = -rho(j)*dv_dx - rho(j)*v(j)*da_dx - v(j)*drho_dx;
%momentum
dvdt_p(j) = -v(j)*dv_dx - (1/gama)*(dTdx+(T(j)/rho(j))*drho_dx);
%energy
dTdt_p(j) = -v(j)*dTdx -(gama-1)*T(j)*(dv_dx+v(j)*da_dx);
%solution updates
v(j)= v(j)+dvdt_p(j)*dt;
T(j)= T(j)+dTdt_p(j)*dt;
rho(j)= rho(j)+drhodt_p(j)*dt;
end
% corrector Method
for j= 2:n-1
dv_dx= (v(j)-v(j-1))/dx;
drho_dx =(rho(j)-rho(j-1))/dx;
dTdx = (T(j)-T(j-1))/dx;
da_dx =(log(a(j))-log(a(j-1)))/dx;
%continuity
drhodt_c(j) = -rho(j)*dv_dx-rho(j)*v(j)*da_dx-v(j)*drho_dx;
%momentum
dvdt_c(j) = -v(j)*dv_dx-(gama^-1)*(dTdx+(T(j)/rho(j))*drho_dx);
%energy
dTdt_c(j) = -v(j)*dTdx -(gama-1)*T(j)*(dv_dx+v(j)*da_dx);
end
%average
dvdt_avg=0.5*(dvdt_c+dvdt_p);
drhodt_avg = 0.5*(drhodt_c+drhodt_p);
dTdt_avg = 0.5*(dTdt_p+dTdt_c);
% Solution update
for i =2:n-1
v(i)= v_old(i)+dvdt_avg(i)*dt;
T(i)= T_old(i)+dTdt_avg(i)*dt;
rho(i)= rho_old(i)+drhodt_avg(i)*dt;
end
%Inlet condition by interpolation
v(1)=2*v(2)-v(3);
%outlet condition by interpolation
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);
%properties
mass = rho.*a.*v;
M = v./((T).^0.5);
p= rho.*T;
%throat values
rho_throat(k)=rho(th);
v_throat(k)=v(th);
T_throat(k)=T(th);
M_throat(k)=M(th);
p_throat(k)=p(th);
mass_throat(k)=mass(th);
%residues at throat.
v_res(k)=dvdt_avg(th);
rho_res(k)=drhodt_avg(th);
T_res(k)=dTdt_avg(th);
%plotting the results
% plotting mass flow rate at different time steps
figure(1)
if k == 50
plot(x, mass, 'b', 'linewidth', 1.5)
hold on
elseif k == 100
plot(x, mass, 'r', 'linewidth', 1.5)
hold on
elseif k == 150
plot(x, mass, 'c', 'linewidth', 1.5)
hold on
elseif k == 200
plot(x, mass, 'g', 'linewidth', 1.5)
hold on
elseif k == 550
plot(x, mass, 'y', 'linewidth', 1.5)
hold on
elseif k == 700
plot(x, mass, 'k', 'linewidth', 1.5)
hold on
elseif k == 1400
plot(x, mass,'m', 'linewidth', 1.5)
grid on
ylabel('mass flow rate');
xlabel('nozzle length (x)');
title('mass flow rate at different time steps Non conservative approach');
legend('50 time steps','100 time steps','150 time steps','200 time steps','550 time steps','700 time steps','1400 time steps');
end
end
figure(3)
plot(linspace(1,nt,nt),v_res);
hold on
plot(linspace(1,nt,nt),rho_res,'g');
hold on
plot(linspace(1,nt,nt),T_res,'k');
grid on
ylabel('residues');
xlabel('timestep');
legend('vplot','rhoplot','Tplot');
title('residual plot');
end
For Conservative Approach:
function[v,rho,T,mass,M,p,v_throat,rho_throat,T_throat,mass_throat,M_throat,p_throat]=conservative_approach(n,L,cfl,nt,th)
% geometry
L=3;
n=31;
x=linspace(0,L,n);
dx= x(2)-x(1);
%Initial Conditions
a=1+2.2*(x-1.5).^2;
th=find(a==1); %thraot node.
rho_initial=ones(1,n);
T_initial=ones(1,n);
for i= 1:n
if (x(i)<=0.5)
rho_initial(i)=1;
T_initial(i)=1;
elseif (x(i)>=0.5 && x(i)<=1.5)
rho_initial(i)=1-0.366*(x(i)-0.5);
T_initial(i)=1-0.167*(x(i)-0.5);
elseif (x(i)>=1.5 && x(i)<=3)
rho_initial(i)=0.634-0.3879*(x(i)-1.5);
T_initial(i)=0.833-0.3507*(x(i)-1.5);
end
end
v = 0.59./(rho_initial.*a);
gamma=1.4;
%Solution term
p= rho_initial.*T_initial;
u1 = rho_initial.*a;
u2 = rho_initial.*a.*v;
u3 = u1.*((T_initial/(gamma-1))+(0.5*gamma*v.^2));
%Time Step calculation based on Courant number.
cfl=0.5;
dt = min(cfl.*(dx./((sqrt(T_initial))+v)));
%Number of iterations
nt=1400;
rho=rho_initial;
T= T_initial;
%time loop
for k=1:nt
%copy of solution terms
u1old=u1;
u2old=u2;
u3old=u3;
%Flux terms
f1= u2;
term=(((gamma-1)/gamma)*(u3-((gamma/2)*((u2.^2)./u1))));
f2=((u2.^2)./u1)+ term;
f3 = ((gamma*u2.*u3)./u1)-(gamma*(gamma-1)*(u2.^3))./(2*u1.^2);
%predictor step
for i =2:n-1
df1dx(i)=(f1(i+1)-f1(i))/dx;
df2dx(i)=(f2(i+1)-f2(i))/dx;
df3dx(i)=(f3(i+1)-f3(i))/dx;
dlogadx(i) = (log(a(i+1))-log(a(i)))/dx;
dadx(i) = (a(i+1)-a(i))/dx;
j2(i)=(1/gamma)*rho(i)*T(i)*dadx(i);
%continuity equation
du1dt_p(i)= -df1dx(i);
%momentum equation
du2dt_p(i)= -df2dx(i)+j2(i);
%energy equation.
du3dt_p(i)= -df3dx(i);
%solution update.
u1(i)=u1(i)+du1dt_p(i)*dt;
u2(i)=u2(i)+du2dt_p(i)*dt;
u3(i)=u3(i)+du3dt_p(i)*dt;
end
%primitive variables
rho = u1./a;
v = u2./u1;
T= (gamma-1)*((u3./u1)-((gamma*v.^2)/2));
%flux terms before corrector step
f1= u2;
term=(((gamma-1)/gamma)*(u3-(gamma*u2.^2)./(2*u1)));
f2=((u2.^2)./u1)+ term;
f3 = ((gamma*u2.*u3)./u1)-((gamma*(gamma-1)/2)*((u2.^3)./(u1.^2)));
%corrector method
for i=2:n-1
df1dx(i)=(f1(i)-f1(i-1))/dx;
df2dx(i)=(f2(i)-f2(i-1))/dx;
df3dx(i)=(f3(i)-f3(i-1))/dx;
dlogadx(i) = (log(a(i+1))-log(a(i)))/dx;
dadx(i) = (a(i)-a(i-1))/dx;
j2(i)=(1/gamma)*rho(i)*T(i)*dadx(i);
%continuity equation
du1dt_c(i)= -df1dx(i);
%momentum equation
du2dt_c(i)= -df2dx(i)+j2(i);
%energy equation
du3dt_c(i)= -df3dx(i);
end
%avergae derivative
du1dt = 0.5*(du1dt_c+du1dt_p);
du2dt = 0.5*(du2dt_c+du2dt_p);
du3dt = 0.5*(du3dt_c+du3dt_p);
%final solution update.
for j =2:n-1
u1(j)=u1old(j)+du1dt(j)*dt;
u2(j)=u2old(j)+du2dt(j)*dt;
u3(j)=u3old(j)+du3dt(j)*dt;
end
%boundary conditions
%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))+(0.5*gamma*v(1)^2));
% outle Boundary
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);
%primitive terms
rho = u1./a;
v = u2./u1;
T= (gamma-1)*((u3./u1)-(gamma/2)*(v.^2));
%properties
mass = rho.*a.*v;
M = v./((T).^0.5);
p= rho.*T;
%throat values
rho_throat(k)=rho(th);
v_throat(k)=v(th);
T_throat(k)=T(th);
M_throat(k)=M(th);
p_throat(k)=p(th);
mass_throat(k)=mass(th);
%residues at throat.
u1_res(k)=du1dt(th);
u2_res(k)=du2dt(th);
u3_res(k)=du3dt(th);
%plotting the results
% plotting mass flow rate at different time steps
figure(2);
if k == 50
plot(x, mass, 'b', 'linewidth', 1.5)
hold on
elseif k == 100
plot(x, mass, 'r', 'linewidth', 1.5)
hold on
elseif k == 150
plot(x, mass, 'c', 'linewidth', 1.5)
hold on
elseif k == 200
plot(x, mass, 'g', 'linewidth', 1.5)
hold on
elseif k == 550
plot(x, mass, 'y', 'linewidth', 1.5)
hold on
elseif k == 700
plot(x, mass, 'k', 'linewidth', 1.5)
hold on
elseif k == 1400
plot(x, mass,'m', 'linewidth', 1.5)
grid on
ylabel('mass flow rate');
xlabel('nozzle length (x)');
title('mass flow rate at different time steps conservative appraoch');
legend('50 time steps','100 time steps','150 time steps','200 time steps','550 time steps','700 time steps','1400 time steps');
end
end
figure(4)
plot(linspace(1,nt,nt),u1_res);
hold on
plot(linspace(1,nt,nt),u2_res,'g');
hold on
plot(linspace(1,nt,nt),u3_res,'k');
grid on
ylabel('residues');
xlabel('timestep');
legend('u1plot','u2plot','u3plot');
title('residual plot');
end
Main Code:
close all
clear all
clc
%input variables
%valve of n can be varied by our own as it is grid size
n = 31;
L = 3;
x = linspace(0,L,n);
%Courent no.
cfl = 0.5;
%nodes at throat
th = ((n-1)/2)+2;
% time steps
nt = 1400;
tic
%Function for non conservative form
[v_nc,rho_nc,T_nc,mass_nc,M_nc,p_nc,v_throat_nc,rho_throat_nc,T_throat_nc,mass_throat_nc,M_throat_nc,p_throat_nc]=non_conservative_approach(n,L,cfl,nt,th);
% Calculating the simulation time
Non_conservative_simulation_time = toc;
tic
% Function for conservative form
[v_c,rho_c,T_c,mass_c,M_c,p_c,v_throat_c,rho_throat_c,T_throat_c,mass_throat_c,M_throat_c,p_throat_c]=conservative_approach(n,L,cfl,nt,th);
%Calculating the simulation time
conservative_simulation_time = toc;
%Analatical Solution Steady State
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);
figure(5)
a1 = subplot(4,1,1);
plot(linspace(1,nt,nt),rho_throat_nc);
hold on
plot(linspace(1,nt,nt),rho_throat_c,'g');
hold on
plot(linspace(1,nt,nt),rho_an,'k');
hold on
grid on
ylabel('Density')
legend('Non-conservative approch','Conseravative approach','Analatical approach')
title('Variation in parameter at throat w.r.t time step')
a2 = subplot(4,1,2);
plot(linspace(1,nt,nt),T_throat_nc);
hold on
plot(linspace(1,nt,nt),T_throat_c);
hold on
plot(linspace(1,nt,nt),T_an);
hold on
grid on
ylabel('Temperature')
legend('Non-conservative approch','Conseravative approach','Analatical approach')
title('Variation in parameter at throat w.r.t time step')
a3 = subplot(4,1,3);
plot(linspace(1,nt,nt),p_throat_nc);
hold on
plot(linspace(1,nt,nt),p_throat_c);
hold on
plot(linspace(1,nt,nt),p_an);
hold on
grid on
ylabel('Pressure')
legend('Non-conservative approch','Conseravative approach','Analatical approach')
title('Variation in parameter at throat w.r.t time step')
a4 = subplot(4,1,4);
plot(linspace(1,nt,nt),M_throat_nc);
hold on
plot(linspace(1,nt,nt),M_throat_c);
hold on
plot(linspace(1,nt,nt),M_an);
hold on
grid on
ylabel('Mach number')
legend('Non-conservative approch','Conseravative approach','Analatical approach')
title('Variation in parameter at throat w.r.t time step')
figure(6)
b1 = subplot(4,1,1);
plot(x,rho_nc);
hold on
plot(x,rho_nc,'*');
hold on
grid on
ylabel('Density')
legend('Non-conservative approch','Conseravative approach','Analatical approach')
title('Variation in properties along X')
b2 = subplot(4,1,2);
plot(x,T_nc);
hold on
plot(x,T_nc,'*');
hold on
grid on
ylabel('Temperature')
legend('Non-conservative approch','Conseravative Solution')
title('Variation in properties along X')
b3 = subplot(4,1,3);
plot(x,p_nc);
hold on
plot(x,p_nc,'*');
hold on
grid on
ylabel('Pressure')
legend('Non-conservative approch','Conseravative Solution')
b4 = subplot(4,1,4);
plot(x,M_nc);
hold on
plot(x,M_nc,'*');
hold on
grid on
ylabel('Mach Number')
legend('Non-conservative approch','Conseravative Solution')
%mass flow rate comparission between two forms of equations
figure(7)
plot(x,mass_nc)
hold on
plot(x,mass_c,'g')
hold on
plot(x,linspace(0.579,0.579,n),'k')
hold on
ylabel('non-dimentionless mass flow rate')
legend('Non-conservative approch','Conseravative approach','Analatical approach')
title('non-dimentionless mass flow rate comparission')
axis([0 3 0.55 0.6])
grid on
y = [rho_nc(th) rho_c(th) 0.634; T_nc(th) T_c(th) 0.833; p_nc(th) p_c(th) 0.528; M_nc(th) M_c(th) 1];
x_axis = categorical({'1.denisty','2.Temperature','3.Pressure','4.mach number'});
figure(10)
labels = {'Non-conservative approch','Conseravative approach','Analatical approach'};
b = bar(x_axis,y,'group');
legend(labels,'location','northoutside','Orientation','Vertical')
a = (1:size(y,1)).';
x_ax = [a-0.25 a a+0.25];
for k = 1:size(y,1)
for m = 1:size(y,2)
text(x_ax(k,m),y(k,m),num2str(y(k,m),'%0.2f'),'HorizontalAlignment','center','VerticalAlignment','bottom')
end
end
Graphs and Explaination
1. Variation of flow properties along the Nozzle length:
Observations from graphs:
From the above graph it is observe that both conservative and non conservative approach are in close agreement at any point in a nozzle.
2. Variation of flow properties with time step at throat
Observations from graph:
a. From the graph we can see that both the conservative and non-conservative tends to steady state after 600th time step.
b. Non-conservative solution is in more close agreement with analytical solution than the conservative solution.
c. The plot that the stability is achieved by the non-conservative form earlier than the conservative form but the amplitude of the initial variations is also higher in the non-conservation form. Also the convergence criteria has been kept equal for both the cases.
3. Mass flow rate variation along X at different time steps-Conservative approach
Observations from graphs
a.It can be seen that the curve for 550,700 and 1400 timesteps are close and lie in the same band
b.The fluctuations are less as we assumed to have constant mass flow rate intially.
4. Mass flow rate variation along X at different time steps non-conservative approach:
Observations from graphs
a.Even here we can see that the curve for 550,700 and 1400 time step lie in the same band are close.
b.The fluctuations are high as the mass flow rate is derived from other parameter.
5. Comparission in mass flow rate:
Observations from graphs
a.Conservative curve is smooth and lies close to the analytical value since mass flow rate is directly calculated and the assumption in the begining is to have to a constant mass flow rate.
b.Non-conservative approaches close to the analytical value only at the throat.
c. the profile of the conservative form is much closer to the exact solution and on the other hand the variation in the mass flow rate of the non-conservative form is larger and also far from the exact analytical solution.
6. Residuals graphs
Observations from graphs
all the resudals from graph for conservative and not conservative is is converging at time step 800 and 600 respectively
Grid Dependency.
The solution is said to be grid indepenent if the solution has no effect even if the number of grid points are increased once a stable solution is oobtained
In our case let us double the no. of grid points to 61 and acoordingly the dt should be changed. This can be obtained from the timestep calculation from courant number for each grid and taking the least.
So for n=61 we should have 2800 time steps, therefore for grid-independency the solution for n=31 must repeat at exactly the double time step.
1. Variation of flow properties along the Nozzle length:
2. Variation of flow properties with time step at throat:
3. Mass flow rate variation along X at different time steps conservative approach:
4. Mass flow rate variation along X at different time steps non-conservative approach:
5. Comparission in mass flow rate
6. Residuals graphs:
Summary:
for n = 31 grid points
for n = 61 grid points
Comparing the plots of grid size n = 31 and n = 61 it can be conclude that the solution is grid independent and solving conservative approach requires more computational effort and thus time than non conservative approach...
CONCLUSION:
conservative non conservative
n = 31 1000 600
n = 61 1300 1000
The flow through the quasi 1-D nozzle is subsonic at the convergent section, transonic at the throat section and supersonic at the divergent section. Among the two forms of the equations, the number of iterations and the time for convergence is lesser in non-conservative form. So it can be concluded that Non-conservative form is faster than conservative form. But conservative form is more accurate than the former in this project.