Simulation of a 1D Super-sonic nozzle flow simulation using MacCormack Method
Introduction:
In this challenge, Navier Stoke's equations along a convergent divergent nozzle is to be solved, and find the variation of different significant parameters. The equations will be solved using MacCormack method, in both conservative form and non conservative form.
As we know that, a flow is possible only if it satisfies the continuity equation, which is derived from law of conservative mass. So, though it says non conservative, it doesn't mean that mass is not conserved. It is only in sense with 'numerical' conservation. It is strictly bound to the numerical techniques. Same holds true for momentum and energy equations.
Conservative form:
The conservative form has flux terms which measures quantities per unit area, if we integrate many such areas, it becomes a volume (control volume), so it basically deals with a fixed control volume. This can also be seen as it contains time rate of change of flux terms which basically means flow between two areas i,e. a fixed control volume. This kind of approach is called as Eulerian approach. It is 'conserving' numerically.
The general form of conservative form of equations is, (or this is how we write to solve it in a conservative form)
`(deltaU)/(deltat) + (deltaF)/(deltax) + (deltaG)/(deltay) + (deltaH)/(deltaz) = J`
Where, U, F, G, H, J are vectors (column matrices). U is called as Solution vector, F, G, H are flux terms or vectors, J is a source term. In case of a transient problem (time marching), the equation is written as,
`(deltaU)/(deltat) = -(deltaF)/(deltax) -(deltaG)/(deltay) -(deltaH)/(deltaz) + J`
Looking at the equation above, we can understand why U is called as a solution vector. It is because it is now a dependent vector (the elements of this column matrix are all dependent variables) and the others are independent vectors. In case of steady state, either F or G or H becomes 'solution vector'. J will be zero in case we neglect body forces.
In our case we want to solve continuity, momentum and energy equations using this method, so we write them in the form as shown above. And solve this by assigning boundary conditions and initial conditions, either explicitly or implicitly (In this case, we are solving it explicitly by using MacCormack method).
After solving for the solution matrix, we will 'decode' it and find the actual parameters, which are called primitive variables.
Now, for the sake of simplicity, using simple mathematics, we will also write flux terms and source terms in terms of elements of solution vector i.e, U(1), U(2) etcetera. This is the basic process to solve a differential equation in Conservative form.
The equations which are going to be solved in this form are,
- Continuity equation:
`(delta(rhoA))/(deltat) + (delta(rhoAV))/(deltax) = 0`
- Momentum equation:
`(delta(rhoAV))/(deltat) + (delta(rhoAV^2 + (1/gamma)pA))/(deltax) = 1/gamma p(deltaA)/(deltax)`
- Energy equation:
`(delta(rho(e/(gamma-1) + gamma/2 V^2)A))/(deltat) + (delta(rho(e/(gamma-1) + gamma/2 V^2)VA + pAV))/(deltax) = 0`
NOTE: The energy here is Internal energy and for an ideal gas, internal energy is a function of temperature (linearly proportional) and can be written as Internal energy (U) = mCdt, where C is specific heat at constant volume, m is mass and dt is a small time step. So we are directly taking temperature in place of internal energy.
The above mentioned equations are now framed into Solution vector, Fllux vector and source term,
U1 = `rhoA`
U2 = `U1*V`
U3 = `U1(e/(gamma-1) + gamma/2 V^2)`
Writing Flux vector in terms of solution vector terms,
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)`
The source is term is also written in terms of elements of solution vector,
J = `(gamma-1)/gamma(U3 - gamma/2(U2^2)/(U1))(delta(ln A))/(deltax)`
NOTE: Writing flux terms and source terms in terms of solution vector terms has no difference theoretically, however on a numerical basis, the big difference is stability.
NOTE 2: All the primitive variable terms here denote non - dimensional form of themselves.
Now, the continuity, momentum and energy equations become,
`(delta(U1))/(deltat) = -(delta(F1))/(deltax)`
`(delta(U2))/(deltat) = -(delta(F2))/(deltax) + J`
`(delta(U3))/(deltat) = -(delta(F3))/(deltax)`
These equations can be solved using MacCormack technique which will be discussed below.
Non Conservative form:
In non conservative form, there are no flux terms like in conservatve form and also we have time rate of change of direct variables, also there are convective terms on the RHS, that too in terms of direct variables, which signifies that it is 'following' the variables across the flow. Which signifies that it has 'moving control volume'. This kind of approach is called as Legrangian approach.
In non conservative form, we directly take the space terms to the right side (dependent variables) and solve for the temporal term (in case of transient state). After assigning initial and boundary conditions, we 'time march' and solve for the dependent term, either explicitly or implicitly. In this case, we are using explicit technique.
The equations in non conservative form are,
- Continuity equation:
`(deltarho)/(deltat) = -rho((deltaV)/(deltax)) - rhoV(delta(lnA))/(deltax) - V (deltarho)/(deltax)`
- Momentum equation:
`(deltaV)/(deltat) = -V (deltaV)/(deltax) - 1/gamma ((deltaT)/(deltax) + T/rho (deltarho)/(deltax))`
- Energy equation:
`(deltaT)/(deltat) = -V(deltaT)/(deltax) - (gamma-1) T ((deltaV)/(deltax) + V (delta(ln A))/(deltax))`
Here too, the primitive varibales are non dimensional forms of themselves.
The two above mentioned forms of equations are solved using MacCormack method, which is discussed below.
MacCormack method:
MacCormack method is a technique to find a time marching solution explicitly. It is a second order accurate in both space and time. The beauty of which is we do not need to solve second order differential equations to get second order accuracy (unlike in Lax - Wendroff method, where we need to derive two times with respect to time). To preserve the second order accuracy, we will have to include both forward and backward differencing schemes in two different steps, which also 'refines' the solution numerically.
We will look at how we proceed in MacCormack method.
Firstly, using Taylor's series expansion, the value of a parameter which is changing with another independent term can be written as,
pi+1=pi+`(deltap)/(deltat) * (trianglet)`
In the above, 'p' is a dependent term which is changing with respect to the independent term 't' (time). The equation physically means that we are adding a 'change' to the previous term for a span of '`trianglet`' (unit time). Now, the differential term in the above term (`(deltap)/(deltat)`) is found out by using numerical discretization techniques, forward and backward, which are first order accurate. And by applying them one after the other, we are maintaining second order accuracy.
Here 't' is time and similarly `trianglet` is a 'time step'.
The above mentioned processes are done in two different steps called as Predictor step and Corrector step.
1. Predictor step:
In Predictor step, we find the time the time rate of change of p by taking all the spatial derivative terms to the RHS and applying forward differencing on them. After solving it we get time rate of change of p corresponding to predictor step.
We will then update the value of p,
pi+1=pi+`(deltap)/(deltat) * (trianglet)` (here, the time rate of change of p is corresonding to predictor step).
2. Corrector step:
on the space terms by having them on the RHS. After solving it, we get time rate of change of p corresponding to corrector step.
We now find the average of time rate of change corresponding to predictor step and also corresponding to corrector step to get the average time rate change.
Average (`(deltap)/(deltat)`) = (Predictor(`(deltap)/(deltat)`) + Corrector(`(deltap)/(deltat)`))/2
Notice that the time rate of change is getting 'refined' by considering values both to the right hand side and to the left hand side.
3. Final solution update:
Now using the average value of time rate change, we find the value of the considered parameter, p
pi+1=pi+`(deltap)/(deltat) * (trianglet)`
NOTE: It is not important to use forward differencing in Predictor step and backward differencing in Corrector step, we can alter them based on the problem.
Time step calculation:
Note that to satisfy stability criteria given by CFL number (Courant Friedrichs Lewy number) (<0.5), we cannot arbitrarily assume a time step value.
The time step value for a non linear partial differential equation is given by,
`trianglet = C (trianglex)/(M + V)`
Here, M is the mach number which is, a non dimensional number, a ratio of speed of an object to the speed of sound. Which is also the ratio of inertial force to elastic force.
The value of M can be taken as `sqrt(Temperature)`
Therefore, time step, `trianglet = C(trianglex)/(sqrtT + V)` as all the terms involved in this are going to be matrices, we can take the minimum value element in the time step matrix.
Sub sonic Super sonic isentropic flow:
The sub sonic super sonic nozzle contains a convergent and divergent section as shown in the figure below,
It is assumed that the nozzle is connected to a reservoir, where the area is infinity (say) and also as the water is still, velocity is also assumed to be zero. At the throat section, where the area is minimum, the mach number is assumed to be one. The mach number, while going upstream, after the throat is less than one (meaning that the speed is less than the speed of sound, hence, subsonic) and going downstream, after the throat is more than one (velocity is greater than the speed of sound, Supersonic). The velocity, pressure and temperature are controlled by the area.
Notice that, according to Bernoulli's equation, total energy should remain constant, which implies that the sum of pressure head, datum head and kinetic head should remain constant and so for this if one head increases, the other decreases. And according to conitnuity equation, which simply is law of conservation of mass, the product `rhoAV` should be constant. Now, in if we look at the nozzle, as the area increases, velocity is also increasing, (if we assume density as constant i.e, the fluid is incompressible i.e, liquids), this is violating the continuity equation and so the flow should not be possible. But if we consider the fluid to be compressible (gases), then density also changes to compensate for the increment in velocity and satisfying the continuity equation and making the flow possible. So for this nozzle, the fluid is a gas.
If we assume the flow to be one dimensional and is changing only in the longitudinal direction and the flow to be isentropic i.e reversible adiabatic, no heat transfer into and out of the system, such type of flow is called as a Quasi static flow.
Boundary conditions:
- From the method of characteristics, at the subsonic entry, two dependent variables need to be specified and the other is to 'float' (Velocity), with respect to time (the velocity at the inside points are calculated by simple linear extrapolation), whereas at the supersonic exit, all variables are allowed to float. All other variables are specified
- In case of supersonic outflow, as all variables float, all the variables are linearly extrapolated
- In conservative and non conservative forms of solving, initial conditions are framed differently, though in conservative method too, we 'fix' two dependent variables at the subsonic entry and let the other float and at the supersonic exit, we let all variables float and find their values by linear extrapolation
Code:
The code to solve the sub sonic super sonic nozzle problem using conservative and non conservative techniques is,
clear all
close all
clc
n = 31;
x = linspace(0,3,n);
dx = x(2) - x(1);
rho = 1-0.3146*x;
T = 1-0.2314*x;
v = (0.1 + 1.09*x).*T.^0.5;
gamma = 1.4;
a = 1+2.2*(x-1.5).^2;
throat = find(a == 1);
m_0_1 = rho.*v.*a;
c = 0.5;
nt = 1000;
for i = 1:n
dt_old(i) = c*(dx/(T(i)^0.5+v(i)));
end
dt = min(dt_old);
[rho1, v1, T1, pressure1, mach1, m1,m_throat1,rho_throat1,v_throat1,pressure_throat1,mach_throat1,T_throat1,m_50_1,m_100_1,m_150_1,m_200_1,m_700_1] = non_conservative(x,dx,rho,a,v,T,dt,nt,gamma,n);
[rho2, v2, T2, pressure2, mach2, m2,m_throat2,rho_throat2,v_throat2,pressure_throat2,mach_throat2,T_throat2,m_50_2,m_100_2,m_150_2,m_200_2,m_700_2,m_0_2] = conservative(x,dx,rho,a,v,T,dt,nt,gamma,n);
%Plots corresponding to non_conservative technique
%Mass flow rate
figure(1)
plot(x,m_0_1,'--k')
hold on
plot(x,m_50_1,'color','r')
hold on
plot(x,m_100_1,'color','b')
hold on
plot(x,m_200_1,'color','m')
hold on
plot(x,m_700_1,'color','g')
xlabel('Length of the nozzle')
ylabel('Mass flow rate')
legend('0dt','50dt','100dt','200dt','700dt')
title('Mass flow rate along the nozzle at different time steps in non conservative technique')
saveas(gcf,'Mass_Flow_rate_diff_NON.png')
%Variation of different properties
figure(2)
subplot(3,2,1)
plot(x,rho1,'r')
ylabel('Density')
xlabel('Length of the nozzle')
subplot(3,2,2)
plot(x,mach1,'m')
ylabel('Mach number')
xlabel('Length of the nozzle')
subplot(3,2,3)
plot(x,v1,'g')
ylabel('Velocity')
xlabel('Length of the nozzle')
subplot(3,2,4)
plot(x,T1,'k')
ylabel('Temperature')
xlabel('Length of the nozzle')
subplot(3,2,5)
plot(x,pressure1,'y')
ylabel('Pressure')
xlabel('Length of the nozzle')
subplot(3,2,6)
plot(x,m1,'b')
ylabel('Mass flow rate')
xlabel('Length of the nozzle')
saveas(gcf, 'Variation_of_diff_prop_non.png')
%Plots corresponding to conservative technique
%Mass flow rate
figure(3)
plot(x,m_0_2,'--k')
hold on
plot(x,m_50_2,'color','r')
hold on
plot(x,m_100_2,'color','b')
hold on
plot(x,m_200_2,'color','m')
hold on
plot(x,m_700_2,'color','g')
xlabel('Length of the nozzle')
ylabel('Mass flow rate')
legend('0dt','50dt','100dt','200dt','700dt')
title('Mass flow rate along the nozzle at different time steps in conservative technique')
saveas(gcf,'Mass_Flow_rate_diff_CON.png')
%Variation of different properties
figure(4)
subplot(3,2,1)
plot(x,rho2,'r')
ylabel('Density')
xlabel('Length of the nozzle')
subplot(3,2,2)
plot(x,mach2,'m')
ylabel('Mach number')
xlabel('Length of the nozzle')
subplot(3,2,3)
plot(x,v2,'g')
ylabel('Velocity')
xlabel('Length of the nozzle')
subplot(3,2,4)
plot(x,T2,'k')
ylabel('Temperature')
xlabel('Length of the nozzle')
subplot(3,2,5)
plot(x,pressure2,'y')
ylabel('Pressure')
xlabel('Length of the nozzle')
subplot(3,2,6)
plot(x,m2,'b')
ylabel('Mass flow rate')
xlabel('Length of the nozzle')
saveas(gcf, 'Variation_of_diff_prop_CON.png')
%Comparison between mass flow rates
figure(5)
plot(x,m1,'color','g')
hold on
plot(x,m2,'color','k')
xlabel('Length of the nozzle')
ylabel('Mass flow rate')
legend('Non Conservative','Conservative')
title('Variation of mass flow rate along the length of the nozzle for non conservative and conservative techniques')
saveas(gcf, 'Variation of mass flow rate.png') Conservative function:
function [rho_new, v, T_new, pressure, mach, m,m_throat,rho_throat,v_throat,pressure_throat,mach_throat,T_throat,m_50,m_100,m_150,m_200,m_700,m_0] = conservative(x,dx,rho,a,v,T,dt,nt,gamma,n)
throat = find(a == 1);
%Initial conditions
for i = 1:n
if(x(i)>=0 && x(i)<=0.5)
rho_new(i) = 1;
T_new(i) = 1;
elseif(x(i)>=0.5 && x(i)<=1.5)
rho_new(i) = 1 - 0.366*(x(i) - 0.5);
T_new(i) = 1 - 0.167*(x(i) - 0.5);
elseif(x(i)>=1.5 && x(i)<=3.5)
rho_new(i) = 0.634 - 0.3879*(x(i) - 1.5);
T_new(i) = 0.833 - 0.3507*(x(i)-1.5);
end
v(i) = 0.59/(rho_new(i)*a(i));
U1(i) = rho_new(i)*a(i);
U2(i) = rho_new(i)*a(i)*v(i);
U3(i) = rho_new(i)*((T_new(i)/(gamma-1))+((gamma/2)*(v(i)^2)))*a(i);
end
m_0 = rho_new.*a.*v;
for k = 1:nt
U1_old = U1;
U2_old = U2;
U3_old = U3;
F1 = U2;
pranav = (((gamma-1)/gamma)*(U3-((gamma/2)*((U2.^2)./U1))));
F2 = ((U2.^2)./U1)+pranav;
F3 = ((gamma*U2.*U3)./U1)-(gamma*(gamma-1)*(U2.^3))./(2*(U1.^2));
%Predictor method
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;
dadx(i) = (a(i+1)-a(i))/dx;
J2(i) = (1/gamma)*(rho_new(i)*T_new(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
%Decoding
for p = 1:n
rho_new(p) = U1(p)/a(p);
v(p) = U2(p)/U1(p);
T_new(p) = (gamma-1)*((U3(p)/U1(p))-((gamma*((v(p))^2))/2));
%Flux terms
F1(p) = U2(p);
pranav(p) = (((gamma-1)/gamma)*(U3(p)-(gamma*(U2(p)^2))/(2*U1(p))));
F2(p) = ((U2(p)^2)/U1(p)) + pranav(p);
F3(p) = (((gamma*U2(p)*U3(p))/U1(p))-(gamma*(gamma-1)*(U2(p)^3))/(2*(U1(p))^2));
%F3(p) = ((gamma*U2(p)*U3(p))/U1(p))-(gamma*(gamma-1)*(U2(p)^3))/(2*U1(p)^2);
end
%Corrected method
for v = 2:n-1
dF1dx(v) = (F1(v)-F1(v-1))/dx;
dF2dx(v) = (F2(v)-F2(v-1))/dx;
dF3dx(v) = (F3(v)-F3(v-1))/dx;
dadx(v) = (a(v)-a(v-1))/dx;
J2(v) = (1/gamma)*(rho_new(v)*T_new(v)*dadx(v));
%Continuity equation
dU1dt_c(v) = -dF1dx(v);
%Momentum equation
dU2dt_c(v) = -dF2dx(v) + J2(v);
%Energy equation
dU3dt_c(v) = -dF3dx(v);
end
%Average values
dU1dt_avg = (dU1dt_c + dU1dt_p)/2;
dU2dt_avg = (dU2dt_c + dU2dt_p)/2;
dU3dt_avg = (dU3dt_c + dU3dt_p)/2;
%Final solution update
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
%Inlet boundaryy conditions
U1(1) = rho_new(1)*a(1);
U2(1) = 2*U2(2)-U2(3);
v(1) = U2(1)/U1(1);
U3(1) = U1(1)*((T_new(1)/(gamma-1))+((gamma/2)*v(1)^2));
%Outlet boundary conditions
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);
%Decoding the primitive terms
rho_new = U1./a;
v = U2./U1;
T_new = (gamma-1)*((U3./U1)-((gamma*v.^2)/2));
%Calculating solutions
m = rho_new.*v.*a;
mach = v./(T_new.^0.5);
pressure = rho_new.*T_new;
%Values at throat
m_throat(k) = m(throat);
rho_throat(k) = rho(throat);
mach_throat(k) = mach(throat);
v_throat(k) = v(throat);
pressure_throat(k) = pressure(throat);
T_throat(k) = T(throat);
%Mass flow rate values at different time steps
if k == 50
m_50 = rho_new.*v.*a;
elseif k == 100
m_100 = rho_new.*v.*a;
elseif k == 150
m_150 = rho_new.*v.*a;
elseif k == 200
m_200 = rho_new.*v.*a;
elseif k == 700
m_700 = rho_new.*v.*a;
end
endNon conservative function:
function [rho, v, T, pressure, mach,m,m_throat,rho_throat,v_throat,pressure_throat,mach_throat,T_throat,m_50,m_100,m_150,m_200,m_700] = non_conservative(x,dx,rho,a,v,T,dt,nt,gamma,n)
throat = find(a == 1);
for k = 1:nt
rho_old = rho;
v_old = v;
T_old = T;
for j = 2:(n-1)
dvdx = (v(j+1) - v(j))/dx;
dlogadx = (log(a(j+1)) - log(a(j)))/dx;
drhodx = (rho(j+1) - rho(j))/dx;
dtdx = (T(j+1) - T(j))/dx;
%Predictor step
%Contuity equation
drhodt_p(j) = -rho(j)*dvdx - rho(j)*v(j)*dlogadx - v(j)*drhodx;
%Momentum equation
dvdt_p(j) = -v(j)*dvdx - (1/gamma)*(dtdx + (T(j)/rho(j))*drhodx);
%Energy equation
dTdt_p(j) = -v(j)*dtdx - (gamma - 1)*T(j)*(dvdx + v(j)*dlogadx);
%Solution update
rho(j) = rho(j) + drhodt_p(j)*dt;
v(j) = v(j) + dvdt_p(j)*dt;
T(j) = T(j) + dTdt_p(j)*dt;
end
for p = 2:n-1
dvdx = (v(p) - v(p-1))/dx;
dlogadx = (log(a(p)) - log(a(p-1)))/dx;
drhodx = (rho(p) - rho(p-1))/dx;
dtdx = (T(p) - T(p-1))/dx;
%Corrector step
%Contuity equation
drhodt_c(p) = -rho(p)*dvdx - rho(p)*v(p)*dlogadx - v(p)*drhodx;
%Momentum equation
dvdt_c(p) = -v(p)*dvdx - (1/gamma)*(dtdx + (T(p)/rho(p))*drhodx);
%Energy equation
dTdt_c(p) = -v(p)*dtdx - (gamma - 1)*T(p)*(dvdx + v(p)*dlogadx);
end
%Finding average values
drhodt_avg = (drhodt_p + drhodt_c)/2;
dvdt_avg = (dvdt_p + dvdt_c)/2;
dTdt_avg = (dTdt_c + dTdt_p)/2;
for l = 2:n-1
%Final solution
rho(l) = rho_old(l) + drhodt_avg(l)*dt;
v(l) = v_old(l) + dvdt_avg(l)*dt;
T(l) = T_old(l) + dTdt_avg(l)*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);
%Calculating solutions
m = rho.*v.*a;
mach = v./((T).^0.5);
pressure = rho.*T;
%Values at throat
m_throat(k) = m(throat);
rho_throat(k) = rho(throat);
pressure_throat(k) = pressure(throat);
T_throat(k) = T(throat);
mach_throat(k) = mach(throat);
v_throat(k) = v(throat);
%Mass flow rate at different time steps
if k == 50
m_50 = rho.*v.*a;
elseif k == 100
m_100 = rho.*v.*a;
elseif k == 150
m_150 = rho.*v.*a;
elseif k == 200
m_200 = rho.*v.*a;
elseif k == 700
m_700 = rho.*v.*a;
end
end
endResults:
- Non conservative form:
Figure 1:
Inferences:
- The density is decreasing along the length of the nozzle. This is to compensate for the increase in velocity for the flow to be possible. We can also feel it physically, as the fluid is becoming lighter (less dense), the kinetic energy is getting increased and so the velocity is increasing
- Velocity is increasing as the internal energy is getting converted to kinetic energy at the throat section, this is because the throat section is kind of an orifice, where the pressure decreases and from Gay Lussac's law for an ideal gas, temperature also decreases. For an ideal gas internal inergy is a direct funtion of temperature, so the internal inergy decreases and because of law of conservation of energy, this energy is getting converted to kinetic energy and so the velocity starts increasing
- The pressure, as mentioned above, is decreasing because of the decrement of area
- Mach number is a direct function of velocity and so as velocity increases, mach number also increases
- Temperature decreases because of decrement in pressure (Gay Lussac's law)
- The mass flow rate is wavy initially and starts decreasing, the least value being at the throat and starts increasing again. This is because in the convergent region, as area and density are decreasing and the velocity is increasing, but the increment in velocity is lower than the decrement in density and area. Also the rate of convergence in convergent portion of the nozzle is more than the rate of increment in the divergent portion of the nozzle (this can be seen in the image of the nozzle). Post the throat region, the increment in velocity and area is more than the decrement in density, so the mass flow rate is decreasing
Figure 2:
Inferences:
- As the flow starts, the mass flow rate is wobbling about the mean steady state (700 time steps), and damps out the disturbance and reaches steady state after 700 time steps. The dashed line is the intial condition (0 time steps)
- As the mass flow rate is the function of area, density and velocity, the disturbances in them (density and velocity) is what's making the wobbling effect. At 700 time steps, the mass flow rate is almost constant, meaning that the density, velocity have reached steady state
- Physically, it can be observed that when the flow starts, a pressure wave is created and gushes out rapidly until it reaches a steady state. The flow is intially turbulent and slowly converts to laminar
- Conservative form:
Figure 1:
Inferences:
- The key difference between conservation and non conservative forms are that in conservative form, flux terms are involved and are changing with respect to the location (convective terms), which is when the volume is fixed (fixed control volume). So, though the equations are same, the basic assumption is different in both conservative and non conservative forms
- The behaviour of the properties will be same as mentioned above except mass flow rate, which is wavy across the length. The disturbances are because of the numerical error of approximation and also such errors which are present in density and velocity, when multiplied, the error is getting amplified
Figure 2:
Inferences:
- It is that initial mass flow rate, at 0 time steps, and when the system reached steady state (700 time steps), is constant and is horizontal. This is because, in case of conservative form, we have flux terms which are measured per unit area. And also flux terms signifies that we are focussed on a particular area and when many such areas are integrated, we get a control volume so basically, conservative form's prime assumption is that the volume is fixed and not moving. This might be the reason for constant mass flow rate initially and also when the system reached steady state
Comparison between mass flow rates in non conservative form and conservative form:
Inferences:
- The mass flow rate in non conservative form is 'wavy'. At the entry and at the exit, the disturbances are more. This might be because the non conservative form inherently assumes that the control volume is not fixed. Whereas, the mass flow rate in conservative form somewhat looks stable across the length which basically assumes that the control volume is fixed (which is best for this case)
- But as time progresses, either of the two methods reaches steady state and the mass flow rate becomes constant then
Conclusion:
- Conservative and non conservative methods both have their own advantages, but whenever the problem we deal with has fixed control volume, it is better to use conservative approach and whenever the problem deals with moving control volume, non conservative approach is better