Week 7 - Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
Problem Statement:
1. To derive the conservation and non-conservation form for the subsonic supersonic flow through the nozzle using Mac-Cormack's technique.
2. To simulate the subsonic and supersonic flow through the nozzle and using both conservative and non-conservative form.
3. To determine the steady-state temperature distribution of flow field variables.
4. To study the difference between the conservative and non-conservative form of governing equations.
Introduction:
In the flow consideration, the nozzle is a convergent-divergent nozzle where the convergent part of the nozzle is attached to the reservoir which has stagnant pressure and temperature. The flow through the nozzle is isentropic flow. The area of the convergent nozzle is large hence the initial velocity is low and at the throat of the nozzle, the flow is sonic where Mach number(M) =1 and with increasing the area the velocity increases at this stage the flow breaks the sound barrier and attains supersonic speed where Mach number>1.
Assumption:
We assume that the flow variables are varying in x-direction ie., along the length of the nozzle. This is done to assume uniform flow properties across the cross-section of the nozzle. This type of flow is called quasi-one-dimensional flow.
Approach:
1. The Governing equations will be written in terms of the partial differential equations suitable for time marching solution for quasi 1D flow.
2. The finite difference model using the Mac Cormack technique is applied to the problem given.
3. Time step size, courant number, and initial and boundary conditions are formulated.
4. Use of the initial conditions to obtain approximate the converge flow variables for the next time step till a steady state is achieved.
Mac Cormack's Method:
In CFD, MacCormack's technique is widely used in the discretization scheme for parabolic partial differential equations. The PDE is of second order.
The MacCormack technique proceeds in two steps i.e., the predictor step followed by the corrector step.
In the predictor step, the gradient value of the governing equation at time step n is estimated with known values at time step n. In this step forward differencing scheme is used.
Then flow variables are calculated at time step n+1 with known values of flow variable at time step n along with gradient.
In the corrector step, the values of the gradient at time step n+1 is estimated with the known values at time step n+1 i .e., obtained from the predicted step. In this step backward differencing scheme is used.
The average of the gradients of the predicted step and corrector step is calculated.
Then flow variables at the n+1 time step are calculated using flow variables at time step n and average gradients are calculated.
Governing Equations:
Continuity Equation:
Momentum Equation:
Energy equation :
Governing Equations in Non-Conservative Form
Continuity Equation:
Momentum Equation:
Energy Equation:
In the above equation (') on the flow variables represent that they are non-dimensionalized.
V' - velocity
T' - temperature
ρ' - is density
A'- is the area of the nozzle for the length x.
Discretization of equations :
We have finally set up that particular form of the equations that will be most appropriate as well as convenient for the time-marching solution of quasi-one-dimensional nozzle flow. This method is a simplified form of solving PDE's, as second-order degree equations are hard to solve,by the use of this method we get second-order accurate time marching without much effort . So, Macormack method is an easy form of explicit finite-difference technique which is second-order-accurate in both space and time .It consists of two steps predictor and corrector.
First, consider the predictor step,In this step, we set up the spatial derivatives as forward differences. Also, to reduce the complexity of the notation, we will drop the use of the prime to denote a dimensionless variable. In what follows, all variables are the nondimensional variables, denoted earlier by the prime notation.As it is a explicit method ,we find the time derivative terms by using the previous state values and update the flow field variables.Then we proceed with corrector step. We then calculate the time derivatives by using these updated values by using backward difference. Then the average value of these time derivatives is found out.
By using this average value we find the corrected values of flow field variables at time t+dt.These steps should be repeated until steady state is reached.
Calculation of time step (Δt) :
Calculate Δt at all the grid points, i = 1 to i = N, by taking C =0.5, and then choose the minimum value among them.The resulting Δt obtained is used in the above equations to calculate flow field variables.
Nozzle Shape and Initial Conditions :
The nozzle shape, A = A(x), is specified and held fixed, independent of time, For the case illustrated in this section, we choose a parabolic area distribution given by
Therefore, in our choice of initial conditions, we have use any knowledge we may have about a given problem in order to intelligently pick some initial conditions so that they are choosen close to steady state values so that there will be no chance of solution to go unstable. For example, in the present problem, we know that density and T decrease and V increases as the flow expands through the nozzle. Hence, we choose initial conditions that qualitatively behave in the same fashion. For simplicity, let us assume linear variations of the flow-field variables, as a function of x. For the present case, we assume the following values at time t = 0.
initial conditions at t=0,
ρ = 1-0.3146*x
T =1-0.2314*x
V = (0.1 + 1.09*x)T^1/2
Boundary conditions :
Another aspect of the numerical solution is that of boundary conditions-an all-important aspect, because without the physically proper implementation of boundary conditions and their numerically proper representation, we have no hope whatsoever in obtaining a proper numerical solution to our flow problem.
Subsonic inflow boundary (point 1) :
Here, we must allow one variable to float.We choose the velocity V1, because on a physical basis we know the mass flow through the nozzle must be allowed to adjust to the proper steady state, and allowing V1 to float makes the most sense as part of this adjustment. The value of V1 changes with time and is calculated from information provided by the flow-field solution over the internal points. We use linear extrapolation from points 2 and 3 to calculate V1 and we take rho(1) =1 and T(1) =1 as fixed values.
V(1) = 2*V(2)-V(3)
Supersonic outflow boundary (point N) :
Here, we must allow all flow-field variables to float We again choose to use linear extrapolation based on the flowfield values at the internal points. Specifically, we have, for the nondimensional variables,
Now we write a program using this data and do simulation to find out the flow field variables and required data.
2) Conservative form:
continuity equation :
Momentum Equation:
Energy Equation:
Now converting these equations to non-dimensional form we get the governing equations as
continuity equation :
Momentum Equation:
Energy Equation:
Let us express the equations for quasi-onedimensional flow in generic form. Let us define the elements of the solutions vector U, the flux vector F, and the source term J as follows.
Now substituting these terms in the above governing equations we get,
To obtain the primitive variables (p, V, T,p, etc.), we must decode the elements U1, U2, and U3 as follows. From the definitions of U1 , U2 , and U3 given above, we have
In the above equations, we note that the flux vector elements F1, F2, and F3 and source terms are couched in terms of the primitive variables. When the computer program is written with F 1 , F2 , and F3 expressed directly in terms of p', V', p', and e', instabilities develop during the course of the time-marching solution .So,we write such that the governing equations are in Pure form in terms of the elements of the solution vector i.e, in terms of U1,U2 and U3.
Finally the next step to do is discretization .We take the governing equations that consists of flux terms and source terms and discretize them by using MacCormack's method.
Predictor step:. To start the calculation, we use the initial conditions for U1, U2, and U3 to calculate the initial values of F1, F2, and F3 and and we use these values to calculate the values of U1,U2 and U3 in the next time step .Calculating the time derivatives by using forward difference to the RHS terms to discretize and these time derivatives are used to calculate the predicted values of U1 ,U2 and U3.
Corrector Step : Now using the predicted values we calculate the new flux terms and source terms .Then we calculate the time derivative terms by using backward difference.Then we calculate the average values of time derivatives of predictor and corrector steps .
Then we calculate the corrected values of U1,U2,U3 by using average values of time derivative.This is done at every node except the boundary nodes since we calculate them by using boundary conditions.Then we also calculate primitive variables.These steps are to be repeated until steady state is obtained.
Initial conditions :
The nozzle shape, A = A(x), is specified and held fixed, independent of time, For the case illustrated in this section, we choose a parabolic area distribution given by
Since the dependent variables being solved in governing are U1, U2, and U3, we need initial conditions for these same variables at time t = 0 in order to start the finite-difference solution. The initial conditions for U1, U2, and U3 also allow initial conditions for F1 F2 and F3 to be obtained from above equations. Such initial conditions for F1, F2, and F3 are needed to form the x derivatives on the right-hand sides of the governing equations at the first time step.
The initial conditions for U1 , U2, and U3 were synthesized by assuming the following variations for ρ' and T':
These variations slightly more realistic than those assumed in the Nozzle Shape and initial conditions for non conservative form this is in anticipation that the stability behaviour of the finite difference formulation using the conservation form of the governing equations might be slightly more sensitive, and therefore it is useful to start with more improve initial conditions than those given for non conservative form. The initial condition for the variation of V' is synthesized by taking advantage of the fact that one of the dependent variables in our governing equations, namely U2 is physically the local mass flow; that is, U2= ρ'A'V'. Therefore,for the initial conditions only, let us assume a constant mass flow through the nozzle and calculate V' as
The value 0.59 is chosen for U2 because it is close to the exact analytical value of the steady-state mass flow (which for this case is 0.579). Therefore, the initial condition for V' as a function of x' is obtained by substituting the p' variation given by above equations.Finally, the initial conditions for U1 , U2, U3 are obtained by substituting the above variations for p',
T', and V' into the definitions above
Boundary conditions :
Subsonic inflow boundary (point 1) :
At the subsonic inflow boundary two properties are held fixed and one is allowed to float, and at the supersonic outflow boundary all properties are allowed to float. In the present formulation, as before, we hold p' and T' fixed at the inflow boundary, both equal to 1.0, and allow V' to float. By holding p' fixed, then U1 at grid point i = 1 is fixed, independent of time, via U1 = p 'A'. That is,
Supersonic outflow boundary (point N) :
At the supersonic outflow boundary all properties are allowed to float.The flow properties at the downstream, supersonic outflow boundary are obtained by linear extrapolation from the two adjacent internal points. If N denotes the grid point at the outflow boundary, then
Now lets write a matlab code which solves the quasi 1D nozzle flow problem using both the conservative and non- conservative forms with the above data.
Matlab code is written to solve the governing equations of Conservative and non-Conservative form using Mac-Cormack's Technique.
Main Script
%Code to Plot the steady state and time wise variation of the flow
%variables of the subsonic supersonic flow through the nozzle and to
%compare the mass flow rate in conservative and non-conservative form.
close all
clear all
clc
%Specifing the number of grid points
n = 31;
%Calling function to calculate the flow variables of the non-conservative form
%of subsonic-supersonic flow through a nozzle.
[mass_nc,rho_nc,t_nc,p_nc,v_nc,a_nc,M_nc] = non_conservative(n);
%Calling for the function to calculate the flow variables of the
%conservative form of subsonic supersonic flow through a nozzle.
[x,mass_c,rho_c,t_c,P_c,v_c,M_c] = conservative(n);
%Plot showing the comparison of mass flow rate of the both forms
figure(6)
plot (x,mass_nc,'g','Linewidth',2)
hold on
plot (x,mass_c,'r','Linewidth',2)
title('Comparison of the mass flow rate of conservative and non conservative form at steady state')
xlabel('length of the nozzle')
ylabel('non dimensional mass flow rate')
legend('Non conservative form','Conservative form')Explanation of the Main script:
1. Defining the number of grid points.
2. Call of the function to calculate the flow variables of the non-conservative form of subsonic supersonic flow through a nozzle that takes an input number of grid points and outputs flow variables.
3. Call of the function to calculate the flow variables of the conservative form of subsonic supersonic flow through a nozzle that takes an input number of grid points and outputs flow variables.
4. Plotting of the graph which compares the mass flow rate of conservative and non-conservative form.
Function to calculate the flow variables of the non-conservative form of subsonic supersonic flow through a nozzle.
% function to calculate the flow variables of the non-conservative form of subsonic supersonic flow through a nozzle.
function [mass,rho,t,P,v,a,M] = non_conservative(n)
%inputs for the non conservation equations
% Total length of Nozzle, L=3
x = linspace(0,3,n);%<-- specifing the values of the position of grid points
dx = x(2)-x(1);%<-- Calculating the grid spacing
%Input Boundary 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;
gamma = 1.4;%<-- heat capacity ratio
throat = find(a==1);%<-- Finding the position of throat where area =1.
%Calculation of time step size.
c=0.5;%<-- CFL criteria
s= sqrt(t);
dt =min(c*dx./(s+v));
%Defining the number of time steps
nt = 1400;
% time loop
for h = 1:nt
%Storing the previous iteration time step converged values.
vold = v;
rhoold=rho;
told = t;
%For loop to calculate the flow variables in the predictor step
for j = 2:n-1
% Inorder to reduce the complexity of the governing equations below
%terms are repeatative and calculated first and then later entered
%into the governing equations.
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;
%Predicted values calculation.
drhodt_p(j) = -rho(j)*dvdx-rho(j)*v(j)*dlogadx-v(j)*drhodx;
dvdt_p(j) = -v(j)*dvdx-(1/gamma)*(dtdx+(t(j)/rho(j))*drhodx);
dtdt_p(j) = -v(j)*dtdx-(gamma-1)*t(j)*(dvdx+v(j)*dlogadx);
%Updating the flow variables with the predicted gradients calculated
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
%Corrector step
%For loop to calculate the flow variables with the predicted flow variables
for jj = 2:n-1
%Corrected values calculation.
dvdx = (v(jj)-v(jj-1))/dx;
dlogadx = (log(a(jj))-log(a(jj-1)))/dx;
drhodx = (rho(jj)-rho(jj-1))/dx;
dtdx = (t(jj)-t(jj-1))/dx;
drhodt_c(jj) = -rho(jj)*dvdx-rho(jj)*v(jj)*dlogadx-v(jj)*drhodx;
dvdt_c(jj) = -v(jj)*dvdx-(1/gamma)*(dtdx+(t(jj)/rho(jj))*drhodx);
dtdt_c(jj) = -v(jj)*dtdx-(gamma-1)*t(jj)*(dvdx+v(jj)*dlogadx);
end
%Taking the average of the gradients of the predicted and corrector
%step
drhodt_avg = 0.5*(drhodt_p+drhodt_c);
dvdt_avg = 0.5*(dvdt_p+dvdt_c);
dtdt_avg = 0.5*(dtdt_p+dtdt_c);
%For loop to calculate the final update solution using the average
%gradients
for ii = 2:n-1
rho(ii) = rhoold(ii) + drhodt_avg(ii)*dt;
v(ii) = vold(ii)+dvdt_avg(ii)*dt;
t(ii) = told(ii)+dtdt_avg(ii)*dt;
end
%Defining initial boundary conditions
v(1) = 2*v(2)-v(3);
%Defining outflow Boundary conditions
rho(n) = 2*rho(n-1)-rho(n-2);
v(n) = 2*v(n-1)-v(n-2);
t(n) = 2*t(n-1)-t(n-2);
%Calculation of flow variable at throat where area(a)=1
rho_t(h) = rho(throat);%<-- density at throat
t_t(h) = t(throat);%<-- Temperature at throat
v_t(h) = v(throat);%<-- velocity at throat
P_t(h) = rho(throat).*t(throat);%<-- Pressure at throat
M_t(h) = v(throat)./sqrt(t(throat));%<-- Mach number at throat
mass = rho.*v.*a;
%if else condition to plot the variation of non dimensional mass flow rate at specific time steps and plot the same
if h == 50
%plotting of mass flow rate at time step 50
figure(7)
plot(x,mass,'Linewidth', 2)
hold on
elseif h == 100
%plotting of mass flow rate at time step 100
plot (x,mass,'Linewidth', 2)
elseif h==150
%plotting of mass flow rate at time step 150
plot(x,mass,'Linewidth', 2)
elseif h==200
%plotting of mass flow rate at time step 200
plot (x,mass,'Linewidth', 2)
elseif h==700
%plotting of mass flow rate at time step 700
plot (x,mass,'Linewidth', 2)
title ('Plot of mass flow rate at different time steps Non conservative form')
xlabel('Nozzle length')
ylabel('Non-dimensional Mass flow rate')
legend ('50','100','150','200','700')
end
end
%Calculation of mach number for the steady state
M = v./sqrt(t);
%Calculation of pressure for the steady state.
P = rho.*t;
%Calculation of mass flow rate for the steady state.
mass = rho.*v.*a;
%Plot for Steady State variation of converged primitive variables inside nozzle
figure(1)
subplot(5,1,1)
plot(x,v,'r','Linewidth',2)
title('Time step:',h)
xlabel('length of nozzle')
ylabel('velcoity')
subplot(5,1,2)
plot(x,t,'y','Linewidth',2)
title('Time step:',h)
xlabel('length of nozzle')
ylabel('Temperature')
subplot(5,1,3)
plot(x,rho,'b','Linewidth',2)
title('Time step:',h)
xlabel('length of nozzle')
ylabel('Density')
subplot(5,1,4)
plot(x,M,'c','Linewidth',2)
title('Time step:',h)
xlabel('length of nozzle')
ylabel('Mach number')
subplot(5,1,5)
plot(x,P,'g','Linewidth',2)
title('Time step:',h)
xlabel('length of nozzle')
ylabel('Pressure')
sgtitle('Steady state of primitive variables for non conservative form')
%time wise variation of flow variables at throat for each time step.
figure(2)
subplot(4,1,1)
plot(t_t)
xlabel('time step')
ylabel('Temperature')
subplot(4,1,2)
plot(rho_t)
xlabel('time step')
ylabel('Density')
subplot(4,1,3)
plot(M_t)
xlabel('time step')
ylabel('Mach number')
subplot(4,1,4)
plot(P_t)
xlabel('time step')
ylabel('Pressure')
sgtitle('Time wise variation of primitive variables for non conservative form')
end Explanation of the above function used for non-conservation form:
1. Function takes the number of grid points(n) as input.
2. With n, the position of the grid points is obtained and grid spacing is calculated.
3. With the position of the grid points are known the initial boundary conditions for the flow variables are calculated.
4. Heat capacity ratio is defined and throat position is found using find command with condition a==1.
5. The governing equation of non-conservative form is of explicit form which is governed by stability criteria where CFL =0.5 with which the time step size is calculated.
6. Time loop is created to calculate the steady-state flow variables. No. of the time loop is of 1400.
7. Inside the time loop, a copy of the previous iteration time step value is made which is later used in the final update of the solution.
8. For the calculation of the predicted values for each time step a for loop is created. The for loop is between the 2nd grid point to the 30th grid point since the initial and last boundary values are defined in boundary conditions.
Flow variables are updated using the predicted step derived above for non-conservative form inside the predictor for loop.
9. Another for loop is created to calculate the gradients for the corrector step using the equation derived for non-conservative form.
10. Then the average of the predictor and corrector gradients are calculated.
11. Another for loop is created for the final update of the flow variables with the gradient average.
12. Inflow and outflow boundary conditions are defined. The reason for defining them at the end of the calculation is that the value of the boundaries depends on the converged values of the inner grid points. Hence the boundary conditions for the flow variables are defined at the end.
13. Storing the values of the flow variables at the throat section to plot the time-wise variation of the flow variables at the throat.
14. If else statements are created to plot the mass flow rate at different time steps to study the variation of the mass flow rate at different time steps till it reaches the steady-state.
15. Mach number(M_nc), pressure(P_nc) and mass flow rate(mass_nc) is calculated.
16. Two plots are created one for steady-state distribution of the flow variables for the length of the nozzle and time-wise variation of the flow variables for the number of time steps.
Function to calculate the flow variables of the conservative form of subsonic supersonic flow through the nozzle.
% function to calculate the flow variables of the
% conservative form of subsonic supersonic flow through a nozzle.
function [x,u2,rho,t,P,v,M] = conservative(n)
%inputs
% Total length of Nozzle, L=3
x = linspace(0,3,n);%<-- specifing the values of the position of grid points
dx=x(2)-x(1);%<-- Calculating the grid spacing
a=1+2.2*(x-1.5).^2;%<-- Defining the area of the convergent divergent nozzle
gamma = 1.4; %<-- heat capacity ratio
%defining the initial conditions
for i = 1:n
if (0<=x(i)&& x(i)<=0.5)
rho(i) =1;
t(i) = 1;
end
if ((0.5<=x(i)&& x(i)<=1.5))
rho(i) = 1-0.366*(x(i)-0.5);
t(i)= 1-0.167*(x(i)-0.5);
end
if (1.5<=x(i)&& x(i)<=3)
rho (i) = 0.634-0.3879*(x(i)-1.5);
t(i) = 0.833 - 0.3507*(x(i)-1.5);
end
end
%Initial conditions
v = 0.59./(rho.*a);
u1 = rho.*a;
u2 = 0.59*ones(1,n);
u3 = u1.*((t./(gamma-1))+(gamma/2)*(v.^2));
throat = find(a==1);%<-- Finding the position of throat where area =1.
%time step size calculations
c=0.5;%<-- CFL criteria
s= sqrt(t);
dt =min(c*dx./(s+v));
nt = 1400;%<-- Defining the number of time steps
% time loop
for h = 1:nt
%Storing the previous iteration time step converged values.
u1old = u1;
u2old = u2;
u3old = u3;
%Calculation of flux variables.
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));
%for loop to calculate the source term j2
for i = 1:n-1
j2(i) = 1/gamma*rho(i)*t(i)*((a(i+1)-a(i))/dx);
end
%Predictor Step
%For loop to calculate the solution variable with flux variable and source term calculated above
for j= 2:n-1
du1dt_p(j) = -(f1(j+1) - f1(j))/dx;
du2dt_p(j) = -(f2(j+1) - f2(j))/dx+j2(j);
du3dt_p(j) = -(f3(j+1) - f3(j))/dx;
%Updating the solution variables with the predicted gradients calculated
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
%Calculation of density and temperature to update the previous iteration values which is used in corrector step
rho = u1./a;
t = (gamma-1)*((u3./u1)-(gamma/2)*(u2./u1).^2);
%Updating the flux variables with the predicted solution variables.
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));
%Corrector step
%For loop to calculate the solution variables with the predicted flux variables
for jj = 2:n-1
du1dt_c(jj) = -(f1(jj) - f1(jj-1))/dx;
du2dt_c(jj) = -(f2(jj) - f2(jj-1))/dx+(1/gamma)*rho(jj)*t(jj)*(a(jj)-a(jj-1))/dx;
du3dt_c(jj) = -(f3(jj) - f3(jj-1))/dx;
end
%Taking the average of the gradients of the predicted and corrector step
du1dt_avg = 0.5*(du1dt_c+du1dt_p);
du2dt_avg = 0.5*(du2dt_c+du2dt_p);
du3dt_avg = 0.5*(du3dt_c+du3dt_p);
%For loop to calculate the final update solution using the average
%gradients
for k = 2:n-1
u1(k) = u1old(k)+ du1dt_avg(k)*dt;
u2(k) = u2old(k)+du2dt_avg(k)*dt;
u3(k) = u3old(k)+du3dt_avg(k)*dt;
end
%Defining initial boundary conditions
u1(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);
%Defining outflow 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);
%Calculation of flow variables at each time step
rho = u1./a;%<-- density calculation
v = u2./u1;%<--velocity calculation
t = (gamma-1)*((u3./u1)-(gamma/2)*v.^2);%<-- temperature calculation
%Calculation of flow variable at throat where area(a)=1
rho_t(h) = rho(throat);%<-- density at throat
t_t(h) = t(throat);%<-- Temperature at throat
v_t(h) = v(throat);%<-- velocity at throat
P_t(h) = rho(throat).*t(throat); %<-- Pressure at throat
M_t(h) = v(throat)./sqrt(t(throat));%<-- Mach number at throat
%if else condition to plot the variation of non dimensional mass flow rate at specific time steps and plot the same
if h == 100
%plotting of mass flow rate at time step 100
figure(5)
plot(x,u2)
hold on
elseif h == 200
%plotting of mass flow rate at time step 200
plot (x,u2)
elseif h==700
%plotting of mass flow rate at time step 700
plot(x,u2)
title ('Plot of mass flow rate at different time steps conservative form')
xlabel('Nozzle length')
ylabel('Non-dimensional Mass flow rate')
legend ('100 time steps','200 time steps','700 time steps')
end
end
%Calculation of converged flow variables for each time step
rho = u1./a;
v = u2./u1;
t = (gamma-1)*((u3./u1)-(gamma/2)*v.^2);
M = v./sqrt(t);
P = rho.*t;
%Plot for Steady State variation of converged primitive variables inside nozzle
figure(3)
subplot(5,1,1)
plot(x,v,'r','Linewidth',2)
title('Time step:',h)
xlabel('length of nozzle')
ylabel('Velocity')
subplot(5,1,2)
plot(x,t,'y','Linewidth',2)
title('Time step:',h)
xlabel('length of nozzle')
ylabel('Temperature')
subplot(5,1,3)
plot(x,rho,'b','Linewidth',2)
title('Time step:',h)
xlabel('length of nozzle')
ylabel('Density ')
subplot(5,1,4)
plot(x,M,'c','Linewidth',2)
title('Time step:',h)
xlabel('length of nozzle')
ylabel('Mach number')
subplot(5,1,5)
plot(x,P,'g','Linewidth',2)
title('Time step:',h)
xlabel('length of nozzle')
ylabel('Pressure')
sgtitle('Steady State of primitive variables for conservative form')
%time wise variation of flow variables at throat for each time step.
figure(4)
subplot(4,1,1)
plot(t_t)
xlabel('time step')
ylabel('Temperature')
subplot(4,1,2)
plot(rho_t)
xlabel('time step')
ylabel('Density')
subplot(4,1,3)
plot(M_t)
xlabel('time step')
ylabel('Mach number')
subplot(4,1,4)
plot(P_t)
xlabel('time step')
ylabel('Pressure')
sgtitle('Time wise variation of primitive variables for conservative form')
end Explanation:
1. Function takes the number of grid points(n) as input.
2. With n, the position of the grid points is obtained and grid spacing is calculated.
3. For loop is created for defining the initial boundary conditions with the equations specified in the introduction for conservative form
4. Since the flow variables are specified in terms of solution variables, the initial boundary conditions are defined for the solution variables and the position of the throat section is also found.
5. The governing equation of non-conservative form is of explicit form which is governed by stability criteria where CFL =0.5 with which the time step size is calculated.
6. Time loop is created to calculate the steady-state flow variables. No. of time loop is of 1400.
7. Inside the time loop, a copy of the previous iteration time step value is made which is later used in the final update of the solution.
8. Flux variables are calculated with solutions variables from each time step.
9. For loop is created in order to calculate the source term. For loop is created because the formula contains the difference term for the area and hence area form 1x30 matrix whereas the density and temperature are 1x31 matrices hence for loop is created which runs for 30 generating 1x30 matrix for j2.
10. For the calculation of the predicted values for each time step a for loop is created. The for loop is between the 2nd grid point to the 30th grid point since the initial and last boundary values are defined in boundary conditions.
Solution variables are updated using the predicted step derived above for the conservative form inside the predictor step for loop.
11. With the predicted solution variable density and temperature and flux variables are calculated.
12. Another for loop is created to calculate the gradients for the corrector step using the equation derived for conservative form with predicted flux variable.
13. Then the average of the predictor and corrector gradients are calculated.
14. Another For loop is created for the final update of the solution variables with the gradient average.
15. Inflow and outflow boundary conditions are defined. The reason for defining them at the end of the calculation is that the value of the boundaries depends on the converged values of the inner grid points. Hence the boundary conditions for the solution variables are defined at the end.
16. Flow variables are found from the converged solution variables for that time step in order to store the flow variables at the throat section for each time loop.
17. If else statements are created to plot the mass flow rate at different time steps to study the variation of the mass flow rate at different time steps till it reaches the steady-state.
18. Final flow variables are calculated at 1400 time steps at the end of the loop.
19. Plot is created for the steady-state and time-wise distribution of the converged flow variables.
Results and observation
Grid independence Test:
The values of flow variables are tabulated at the throat section of the nozzle.
| Flow Variables | Analytical solution | 31 Grid points | 61 grid points | ||
| Conservation form | Non-Conservation form | Conservation form | Non-Conservation form | ||
| Density (ρ">ρρ) | 0.634 | 0.652 | 0.638 | 0.639 | 0.637 |
| Temperature(T) | 0.833 | 0.841 | 0.836 | 0.835 | 0.835 |
| Pressure(P) | 0.528 | 0.548 | 0.534 | 0.534 | 0.533 |
| Mach Number(M) | 1 | 0.978 | 0.999 | 0.998 | 0.997 |
From the above table, we infer that the values of the flow variables changes in the order of `10^-2` for 61 number of grid points when compared to 31 number of grid points. The non-Conservation form results in solutions closer to the analytical solution this is because the rate of decay of residuals for the non-conservation form is higher than the conservation form hence results in a slightly accurate solution than the conservation form
Output Plots
Non-Conservation form
- Steady-state distribution of primitive variables with respect to the length of the nozzle.
- Time-wise distribution of flow variables at different time steps\
- Mass Flow rate at different time steps
Conservation Form
- Steady State distribution of flow variables along the nozzle length
Time-wise variation of flow variables at different time steps
Mass flow rate at different time steps:
Mass flow comparison of Conservation and non Conservation form at 1400 time steps
Explanation of the Output graph:
The flow variables density, temperature, and pressure decrease along the length of the nozzle whereas velocity and Mach number increases along the length of the nozzle, this type of behavior of the flow are because in the converging part of the nozzle the Mach number is less than 1 as it approaches throat, the area of the nozzle decreases and at the throat, the Mach number is equal to 1. The flow expands isentropically, as a result, the flow accelerates in the divergent part of the nozzle also due to an increase in area and decrease in pressure which results in a decrease in density.
Variations in flow variables are high at the initial time steps because the time derivatives are large and oscillate. The oscillations are due to unsteady compression and expansion of waves which propagate through the nozzle. As the time steps increase the flow reaches a steady-state the time derivatives become small and change in the order of six in their magnitude at larger time steps. As the flow attains the steady-state the time derivative becomes zero but when numerically solved this is not achieved. These non zero values are called residuals which are used by CFD experts in comparing the time marching solution of the steady-state with their magnitude and the rate of decay.
The same effect can be used to explain the non-dimensional mass flow rate at initial times there are large oscillations of mass flow rate along the length of the nozzle as the flow approaches steady state the mass flow rate is constant which can be seen at a 700th-time step it approaches the steady-state value 0.58 which is close to analytical value.
In the comparison of the mass flow rate between conservation and non-conservation form plot, we can infer that the mass flow rate is closer to the analytical solution. With this, we can say that the conservation form of the governing equations performs well in conserving the mass flow rate than the non-conservation form this is because in the conservation form the mass flow rate is one of the dependent variables whereas in the non-conservation form mass flow rate is calculated from the dependent variables i.e., flow variables.
Conclusion:
1. In conservation form mass flow rate is more conserved which yields a better mass flow distribution whereas non-conservation form results in better-converged solution closer to the analytical solutions this is because of greater decay of the residuals.
2. Non-Conservation form results in accurate flow variables whereas the Conservation form results in accurate flux variables.
3. Computational effort needed for conservation form is more as the flow variables need to be calculated from the solution variables whereas in the non-conservation form the converged flow variables are calculated directly.
4. Conservation form of governing equations is used in shock capturing whereas non-conservation form cannot be applied.
Shock is a sudden change in the flow variables in the direction of flow. The conservation form of governing equations is more suitable in fluid flows that have shocks involved because the primary dependent variables in the equation are the solution variables which results in smooth and stable results whereas in non-conservation form the primary dependent variables are flow variables which result in discontinuous solution during shocks. Hence Conservation form is better than the non-conservation form which involves shock-capturing.
From the above discussion, we cannot conclude the conservative form of governing equations is superior to the non-conservative form or vise versa. The type of governing equations is chosen based on the problem at hand.