Simulation of a 1D Super-sonic nozzle flow simulation using MacCormack Method
Aim:
To simulate the isentropic quasi 1-D nozzle flow using the MacCormack's Method for both conservative and non-conservative forms of the governing equations and compare the results regarding the primitive variables such as Temperature, Pressure, Density, and Mass Flow rate.
Theory:
Consider a Convergent-Divergent nozzle as the one shown in the figure below.
We will consider the flow through the nozzle here as isentropic. The flow comes from a reservoir where the Area of the reservoir is very large ( ideally `A rightarrow oo`) and hence the velocity is very small. The flow eventually expands to supersonic speeds beyond the throat achieving the Mach-Number greater than 1. As we will see the local flow around the throat is subsonic tending to supersonic i.e. the Mach Number is almost 1. We denote our constant values in the throat section as `v_0 , T_0 and P_0`. To solve the flow field equations the MacCormack's Method is used.
Governing Equations
We derive the governing equations for the flow field using three basic principles in Mechanics:
- Mass Conservation
- Newton's Second Law of Motion
- Energy Conservation
Consider a Control Volume as sketched in the figure shown below :
As the control area is itself a function of the x, we must start the derivation with the 1-D Integral forms and work our way through there. Consider the enlarged view of the control area with the flow field variables as shown:
Continuity Equation:
Consider the 1-D Integral form of the continuity equation:
`int int int_V frac(partialrho)(partialt)dV + int int_S rho v*dS = 0`
where V denotes the control volume and S denotes the control surface area. Expanding the integral we get:
`frac(partial)(partialt)intintint_V rho dV + int int_Srhov.dS = 0`
Now, we get
`frac(partial)(partialt)int int_VrhodV = frac(partial)(partialt)(rho A dx)`
and,
`int int_S rhov.dS = (rho+drho)(v+dv)(A+dA) - rhovA`
Exoanding the terms on the right we get:
`int int_Srhov.dS = rho v A +rho v dA +rhodvA + vAdrho + vdrhodA + Adrhodv + rhodAdv - rhoAv`
Simplifying the terms and neglecting the higher order differential terms with multiplication of two differentials we get,
`int int_Srhov.dS = rhoAdv + Avdrho+rhovdA = d(rhovA)`
Now, combining the terms we get
`frac(partial)(partialt)(rhoAdx) + d(rhoAv) = 0`
Dividing throughout by dx we get,
`frac(partial)(partialt)(rhoA) + frac(partial)(partialx)(rhoAv) = 0`
which is the conservative form of the continuity equation.
Momentum Equation:
Consider the integral form of the momentum equation as shown below:
`frac(partial)(partialt)int int int_V (rhou)dV + int int_S(rhouv).dS = -int int_S(P.dS)_x`
This equation involves the calculation of three terms. We shall do the calculation individually for the three terms and will combine them:
`frac(partial)(partialt)int int int_V(rhou)dV = frac(partial)(partialt)(rhovAdx)`
Now consider the second term:
`int int_S(rhouv).dS = (rho+drho)(v + dv)^2(A+dA)-(rhov^2A)`
Simplifying the terms and neglecting the higher order differential terms we get,
`int int_S(rhouv).dS = v^2Adrho + 2rhoAvdv + rhov^2dA + (H.O.T)`
which can be simplified further to:
`int int_S(rhouv).dS = d(rhoAv^2)`
Now consider the third term:
`int int_s(P.dS)_x`. Consider the following diagram to understand the terminology properly:
As can be seen, the y-components across the slant faces will get canceled out and we will get:
`int int_s(P.dS)_x = (p+dP)(A+dA) - PA - 2Pfrac(dA)(2)`
Simplifying the terms we get:
`int int_s(P.dS)_x = AdP`
Combining the terms we get:
`frac(partial)(partialt)(rhovA) + frac(partial)(partialx)(rhov^2A) = -Afrac(partialP)(partialx)`
This is the conservative form of the momentum equation:
Energy Equation:
Consider the integral form of the energy equation as shown below:
`frac(partial)(partialt)intintint_Vrho(e + frac(v^2)(2))dV + intint_Srho(e + frac(v^2)(2))v.dS = -intint_S(Pv).dS`
Again to simplify this equation we shall look at each term separately:
term1: `frac(partial)(partialt)intintint_Vrho(e + frac(v^2)(2))dV`
`frac(partial)(partialt)intintint_Vrho(e + frac(v^2)(2))dV = frac(partial)(partialt)(rho[e + frac(v^2)(2)]Adx)`
term2: `intint_Srho(e + frac(v^2)(2))v.dS`
`intint_Srho(e + frac(v^2)(2))v.dS = (rho+drho)(e +de+ frac((v+dv)^2)(2))(v+dv)(A+dA) - rho(e + frac(v^2)(2))vA`
Simplifying the terms and neglecting the higher order differential terms we get,
`intint_Srho(e + frac(v^2)(2))v.dS = frac(d(rhoAv^3))(2) + d(rhoAve)`
term3:`intint_S(Pv).dS`
`intint_S(Pv).dS = (P+dP)(v+dv)(A+dA) - PvA - 2Pvfrac(dA)(2)`
Simplifying the terms and neglecting the higher order differential terms we get,
`intint_S(Pv).dS = d(PAv)`
Combining the terms we get:
`frac(partial)(partialt)(rho[e + frac(v^2)(2)]Adx) + frac(d(rhoAv^3))(2) + d(rhoAve) = -d(PAv)`
`frac(partial)(partialt)(rho[e + frac(v^2)(2)]Adx) + d(rhoAv[e + frac(v^2)(2)])= -d(PAv)`
Dividing throughout by dx we get:
`frac(partial)(partialt)(rho[e + frac(v^2)(2)]A) + frac(partial)(partialx)(rhoAv[e + frac(v^2)(2)])= -frac(partial)(partialx)(PAv)`
This is the conservative form of the momentum equation:
Now to solve these equations numerically we will need a technique. Here we shall be implementing the MacCormack's Method to solve these flow field equations numerically.
MacCormack's Method
MacCormack's method is a variant of the Lex-Wendroff approach but is much simpler than that. The basic idea in the MacCormack's method is to predict an intermediate value of the differential term and then using that correcting the terms and then adding it to the value from the previous time step so as to get the value at the next time step. For illustration purposes consider the simple equation:
`frac(partial)(partialt)B + frac( partial)(partialx)B = 0`
where T and B are some variables we want to solve for. It can be written as:
`frac(partial)(partialt)B = - frac( partial)(partialx)B`
Now in the predictor method we forward difference the spatial terms as shown below:
`(frac(partialB)(partialt))_p = -(frac(B_(i+1) - B_i)(Deltax))`
where the right-hand terms are known from the previous time step. Now after performing the predictor step it is time we perform the corrector step. Note that in corrector step we must use the values as obtained from predictor step i.e.
`\bar(B) = B + (frac(partialB)(partialt))_p dt`
We will use the value of `bar(B)` to perform the corrector step.
It is imperative to note that in the corrector step we use the backward differencing for the spatial terms as shown below:
`(frac(partialB)(partialt))_c = -(frac(B_(i) - B_(i-1))(Deltax))`
Then we compute the average time derivative as :
`(frac(partialB)(partialt))_(av) = frac((frac(partialB)(partialt))_(p)+ (frac(partialB)(partialt))_(c))(2)`
And then we compute the next time step value as :
`B^(n+1) = B^n +(frac(partialB)(partialt))_(av) dt`
After the cmpletion of this step one iteration of the MacCormack Method is performed. This process is continued with `B^(n+1)` as the value for the next time step value calculation untill the desired accuracy or the time is concerned.
Now with enough theory, we are in a position to tackle the problem:
CONSERVATIVE APPROACH
The conservative form of the equations are :
`frac(partial)(partialt)(rhoA) + frac(partial)(partialx)(rhoAv) = 0`
`frac(partial)(partialt)(rhovA) + frac(partial)(partialx)(rhov^2A) = -Afrac(partialP)(partialx)`
`frac(partial)(partialt)(rho[e + frac(v^2)(2)]A) + frac(partial)(partialx)(rhoAv[e + frac(v^2)(2)])= -frac(partial)(partialx)(PAv)`
Now to start with the MaCormack Method we need to seperate the time derivative terms from the spatial derivative terms:
`frac(partial)(partialt)(rhoA) =- frac(partial)(partialx)(rhoAv)`
`frac(partial)(partialt)(rhovA) =- frac(partial)(partialx)(rhov^2A) -Afrac(partialP)(partialx)`
`frac(partial)(partialt)(rho[e + frac(v^2)(2)]A) = - frac(partial)(partialx)(rhoAv[e + frac(v^2)(2)]) -frac(partial)(partialx)(PAv)`
We can take the energy as `e = C_vT` , where T is the temperature. We can also write the Pressure as `P = rhoRT`. Substituting the values and somplifying we get:
`frac(partial)(partialt)(rhovA) =- frac(partial)(partialx)(rhov^2A) -Afrac(partial(rhoRT))(partialx)`
`frac(partial)(partialt)(rho[C_vT + frac(v^2)(2)]A) = - frac(partial)(partialx)(rhoAv[C_vT + frac(v^2)(2)]) -frac(partial)(partialx)(rhoRTAv)`
Now we need to non-dimensionalize the equations in order to facilitate easy solving
We choose the non-dimensional parameters as :
- `x' = frac(x)(L)` where L is the length of the nozzle
- `t' = frac(t)(L//v_0)` where `v_0` is the sonic speed
- `rho' = rho/rho_0`
- `T' = T/T_0`
- `v' = v/v_0` and finally
- `A' = A/A_0`
where the values `rho_0, T_0,v_0 and A_0` are the values at the throat of the nozzle.
Substituting these variables in the continuity equation we get:
`frac(partial)(partial(t'L//v_0))(rho'rho_0A'A_0) =- frac(partial)(partial(x'L))(rho'rho_0A'A_0v'v_0)`
`frac(partial)(partial(t'))(rho'A') =- frac(partial)(partial(x'))(rho'A'v')`
Substituting these variables now into the momentum equation we get:
`frac(partial)(partial(t'L//v_0))(rho'rho_0v'v_0A'A_0) =- frac(partial)(partial(x'L))(rho'rho_0(v')^2v_0^2A'A_0) -A'A_0frac(partialrho'rho_0RT'T_0)(partial(x'L))`
Using the fact that `v_0^2 = gammaRT_0`, we can simplify the equation as :
`rho_0v_0^2A_0/Lfrac(partial)(partial(t'))(rho'v'A') =- rho_0v_0^2A_0/Lfrac(partial)(partial(x'))(rho'(v')^2A') -A_0rho_0RT_o/LA'frac(partial(rho'T'))(partial(x'))`
`frac(partial)(partial(t'))(rho'v'A') =- frac(partial)(partial(x'))(rho'(v')^2A') -1/gammaA'frac(partial(P'))(partial(x'))`
Making a substitution as :`1/gammafrac(partial(P'A'))(partialx') = 1/gammaP'frac(partial(A'))(partialx') + 1/gammaA'frac(partial(P'))(partialx')`
We can write the final equation as :
`frac(partial)(partial(t'))(rho'v'A') =- frac(partial)(partial(x'))(rho'(v')^2A' ) -1/gammaP'frac(partial(A'))(partialx') + 1/gammafrac(partial(P'A'))(partialx')`
Simplifying it furthur we get :
`frac(partial)(partial(t'))(rho'v'A') =- frac(partial)(partial(x'))(rho'(v')^2A' + 1/gammaP'A' ) +1/gammaP'frac(partial(A'))(partialx')`
Substituting these variables into the energy equation we get:
Energy Equation: `frac(partial)(partialt)(rho[C_vT + frac(v^2)(2)]A) = - frac(partial)(partialx)(rhoAv[C_vT + frac(v^2)(2)]) -frac(partial)(partialx)(PAv)`
Non-Dimensionalising the equation we get:
`frac(partial)(partial(t'L//v_0))(rho'rho_0[C_vT'T_0 + frac((v')^2v_0^2)(2)]A'A_0) = - frac(partial)(partial(x'L))(rho'rho_0A'A_0v'v_0[C_vT'T_0 + frac((v')^2v_0^2)(2)]) -frac(partial)(partial(x'L))(P'P_0A'A_0v'v_0)`
We also know that `C_v = R/(gamma-1)`. Using this relation we get:
`frac(partial)(partialt')(rho'[R/(gamma-1)T' + frac((v')^2gamma)(2)]A') = - frac(partial)(partial(x'))(rho'A'v'[R/(gamma-1)T' + frac((v')^2gamma)(2)]) -frac(partial)(partial(x'))(P'A'v')`
Now after the non-dimensionalization of the equations we have three equations that we want to solve :
- `frac(partial)(partial(t'))(rho'A') =- frac(partial)(partial(x'))(rho'A'v')`
- `frac(partial)(partial(t'))(rho'v'A') =- frac(partial)(partial(x'))(rho'(v')^2A' + 1/gammaP'A' ) +1/gammaP'frac(partial(A'))(partialx')`
-
`frac(partial)(partialt')(rho'[R/(gamma-1)T' + frac((v')^2gamma)(2)]A') = - frac(partial)(partial(x'))(rho'A'v'[R/(gamma-1)T' + frac((v')^2gamma)(2)]) -frac(partial)(partial(x'))(P'A'v')`
Now as it can be seen in the non-dimensional equations the primitive variables are coupled and hence we must devise a new strategy to discretize the equations. We define the following parameters :
The state variables:
- `U_1 = rho'A'`
- `U_2 = rho'A'v'`
-
`U_3 = rho'((T')/(gamma-1) + gamma/2(v')^2)A'`
The flux variables:
- `F_1 = rho'A'v'`
- `F_2 = rho'A'(v')^2 + 1/gammaP'A'`
- `F_3 = rho'((T')/(gamma-1) + gamma/2(v')^2)A'v' + P'A'v'` and
- `T_2 = 1/gammaP'frac(partialA')(partialx')`
From here onwards we shall denote the non-dimensional variables as `P,v,T,rho` etc. for simplicity reasons
Using all the defined terms we can write the equations as:
`frac(partialU_1)(partialt) = -frac(partialF_1)(partialx)`
`frac(partialU_2)(partialt) = -frac(partialF_2)(partialx) + J_2`
`frac(partialU_3)(partialt) = -frac(partialF_3)(partialx)`
It can be shown very easily that the Flux terms are related to the state variables as :
`F_1 = U_2`
`F_2 = U_2^2/U_1 + (gamma-1)/gamma(U_3 - gamma/2U_2^2/U_1)` and
`F_3 = U_2U_3/U_1 + U_2(gamma-1)[U_3/U_1 - gamma/2U_2^2/U_1^2]`
These equations will be used to calculate the state variables in between the predictor and corrector step and after one iteration sweep.
Now we can apply the MacCormack's method to these equations :
Predictor Step:
`(frac(partialU_1)(partialt))_p = -(frac(F_1(i+1)-F_1(i))(Deltax))`
`(frac(partialU_2)(partialt))_p = -(frac(F_2(i+1) - F_2(i))(Deltax)) + J_2(i)`
`(frac(partialU_3)(partialt))_p = -(frac(F_3(i+1)-F_3(i))(Deltax))`
Now we can update the intermediate values as :
`bar(U_1) = U_1 + (frac(partialU_1)(partialt))_p`
`bar(U_2) = U_2 + (frac(partialU_2)(partialt))_p` and
`bar(U_3) = U_3 + (frac(partialU_3)(partialt))_p`
Using these values the intermediate Flux Values can be calculated as discussed above:
So now `bar(F_1) , bar(F_2) and bar(F_3) and bar(J_2)` are knwon.
Now using these intermediate values we apply the corrector step as :
`(frac(partialU_1)(partialt))_c = -(frac(\bar(F_1)(i+1)-(\barF_1)(i))(Deltax))`
`(frac(partialU_2)(partialt))_c = -(frac(\barF_2(i+1) - \barF_2(i))(Deltax)) + \barJ_2(i)`
`(frac(partialU_3)(partialt))_c = -(frac(\barF_3(i+1)-\barF_3(i))(Deltax))`
After the completion of the corrector step we move to the final solution Update:
`U_1^(n+1) = U_1^n + frac( (frac(partialU_1)(partialt))_c + (frac(partialU_1)(partialt))_p)(2)`
`U_2^(n+1) = U_2^n + frac( (frac(partialU_2)(partialt))_c + (frac(partialU_2)(partialt))_p)(2)`
`U_3^(n+1) = U_3^n + frac( (frac(partialU_3)(partialt))_c + (frac(partialU_3)(partialt))_p)(2)`
After the completion the final primitive variables are calculated as discussed above:
Boundary Conditions:
For the Boundary Conditions `U_1` remains unchanged ,
`U_2(i=1) = 2U_2(i=2) - U_2(i=3)`
and similarly `U_3` is calculated
For the end point boundary conditions:
`U_1(i=n) = 2U_1(i=(n-1)) - U_1(i=(n-2))`
`U_2(i=n) = 2U_2(i=(n-1)) - U_2(i=(n-2))`
`U_3(i=n) = 2U_3(i=(n-1)) - U_3(i=(n-2))`
Assumption of the initial Profiles:
Assume :
`rho = {(1 | 0\leqxleq0.5),(1-0.366(x-0.5) | 0.5\leqxleq1.5),(0.634-0.3879(x-1.5) |1.5\leqxleq3):}`
and
`T = {(1 | 0\leqxleq0.5),(1-0.167(x-0.5) | 0.5\leqxleq1.5),(0.833-0.3507(x-1.5) |1.5\leqxleq3):}`
The initial Area profile is assumed as :
`A = 1 + 2.2(x-1.5)^2 | 0\leqx\leq3`
The following MATLAB CODE is implemented for initialising the variables to be solved using the MacCormack method:
clear;
close all;
clc;
%
% In this program I shall be coding the nozzle problem in 1-D using the
% MacCormack's method with the use of the predictor corrector method
%
% The problem domain is selected as [0,3] for the non-dimensional variable
% x
%
% The number of grid points used for the discretization including the
% boundary points are n = 31
%
% The initial temperature profiles are assumed @t=0
%
% The equations gorverning the nozzle flow are non-linear PDE and are
% conservative in nature
%
% Defining the number of grid points
n = 31;
% Defining the problem domain
xx = linspace(0,3,n);
% Defining the grid size
dx = xx(2) - xx(1);
% Define the gamma
gamma = 1.4;
% Defining the Area profile
A = 1 + 2.2*(xx-1.5).^2;
% Defining the peicewise continuous profile for density
for j=1:length(xx)
% for x less than 0.5
if xx(j)<=0.5
% Defining the profile
rho_initial(j) = 1;
% for x in between 0.5 and 1.5
elseif xx(j)<=1.5
% Defining the profile
rho_initial(j) = 1-0.366*(xx(j)-0.5);
% for the remaining domain
else
% Define the profile
rho_initial(j) = 0.634 - 0.3879*(xx(j)-1.5);
end
end
% Defining the peicewise continuous profile for Temperature
for j=1:length(xx)
% for x less than 0.5
if xx(j)<=0.5
% Defining the profile
T_initial(j) = 1;
% for x in between 0.5 and 1.5
elseif xx(j)<=1.5
% Defining the profile
T_initial(j) = 1-0.167*(xx(j)-0.5);
% for the remaining domain
else
% Define the profile
T_initial(j) = 0.833 - 0.3507*(xx(j)-1.5);
end
end
% Defining the profile of velocity from the U2
v_initial = 0.59./(rho_initial.*A);
% Defining the tolerance
tol = 1e-6;
% Defining the initial Pressure
P_initial = rho_initial.*T_initial;
% Defining the time interval for movie plot
interval_movie_plot = 100;
% Defining the point at which we want to study the variation
point = round(n/6);
Explanation: This function basically generates the initial variables to be used in solving the nozzle flow problem like the initial profiles of the primitive variables, the number of grid points to be used, the interval length of plotting the movie and so on.
We shall see the initial plot of the state variables which is implemented using the following MATLAB Code:
%% Plotting the initial conditions
figure(1);clf;
set(gcf,'Position',[0 0 1600 1600]);
subplot(221);
plot(xx,rho_initial,'r-','linew',2);
xlabel('Non-Dimensional X-Axis');
ylabel({'Non-Dimensional';' Density'});
title('Initial Density');
set(gca,'fontsize',17);
grid on;
grid minor;
subplot(222);
plot(xx,v_initial,'r-','linew',2);
xlabel('Non-Dimensional X-Axis');
ylabel({'Non-Dimensional';' Velocity'});
title('Initial Velocity');
set(gca,'fontsize',17);
grid on;
grid minor;
subplot(223);
plot(xx,T_initial,'r-','linew',2);
xlabel('Non-Dimensional X-Axis');
ylabel({'Non-Dimensional';' Temperature'});
title('Initial Temperature');
set(gca,'fontsize',17);
grid on;
grid minor;
subplot(224);
plot(xx,P_initial,'r-','linew',2);
xlabel('Non-Dimensional X-Axis');
ylabel({'Non-Dimensional';' Pressure'});
title('Initial Pressure');
set(gca,'fontsize',17);
grid on;
grid minor;
The following Initial plots are obtained :
A solver is built in MATLAB which takes the initial condition from the above code and solves the quasi 1-D nozzle flow problem. The code implemented to construct the solver is depicted here:
function [rho_final,v_final,T_final,P_final,m_rate,r,v,T,P,m,...
iterations_to_steady,M_1,M_2,Time_taken] = ...
conservative_MacCormack_solver(gamma,...
rho_initial,T_initial,A,P_initial,tol,xx,n,a,point)
%
% This function will take input as the initial density and initial
% Temperature as the input and the steady state tolerance desired and will
% produce M_1 and M_2 movie loop counters holding the frames, the primitive
% state variables at the final step and the variation of non-dimensional
% variables with iterations at the grid point specified by the user
%
% The conservative equations of the 1-D quasi flow are given by:
%
% d/dt(U1) = - d/dx(F1)
% d/dt(U2) = - d/dx(F2) + J2
% d/dt(U3) = - d/dx(F3)
%
% Defining the grid spacing
dx = xx(2)-xx(1);
% Defining the time step size
dt = 0.0267;
% Time derivative terms
U1_initial = rho_initial.*A;
U2_initial = 0.59*ones(numel(xx),1)';
% Defining the profile of velocity from the U2
v_initial = U2_initial./(rho_initial.*A);
U3_initial = rho_initial.*( T_initial/(gamma-1) + ...
(gamma/2)*v_initial.^2 ).*A;
%
% The spatial terms
%
F1_initial = U2_initial;
F2_initial = rho_initial.*A.*v_initial.^2 + ...
(1/gamma)*P_initial.*A;
F3_initial = rho_initial.*( T_initial/(gamma-1) + ...
(gamma/2).*v_initial.^2 ).*A.*v_initial + ...
P_initial.*A.*v_initial;
% Defining the movie loop counter
mm = 0;
% Defining the state error
state_error = 9e9;
% Initialising the loop time counter
nt = 0;
%
% It is imperative to note that these variables are non-dimensional
%
% Now we begin our MacCormack method
%
% Main outer time loop
while(state_error>tol)
% Increement the loop counter
nt = nt + 1;
% Keep a copy of the old variables
rho_old = rho_initial;
v_old = v_initial;
T_old = T_initial;
P_old = P_initial;
U1_old = U1_initial;
U2_old = U2_initial;
U3_old = U3_initial;
F1_old = F1_initial;
F2_old = F2_initial;
F3_old = F3_initial;
% Predictor Method
%
% First we shall calculate the predictor variables using the forward
% differencing scheme
for j = 2:n-1
% Calculating the individual terms in the equations so as to easy
% to program and easy to debug
%
% The spactial terms using the forward differencing
dF1_dx = ( F1_initial(j+1) - F1_initial(j) )/dx;
dF2_dx = ( F2_initial(j+1) - F2_initial(j) )/dx;
dF3_dx = ( F3_initial(j+1) - F3_initial(j) )/dx;
dlogA_dx = ( log(A(j+1)) - log(A(j)) )/dx;
J2_initial = (gamma-1)/(gamma)*( U3_initial(j) - (gamma/2) *...
U2_initial(j)^2/U1_initial(j) ) * dlogA_dx ;
% Performing the predictor method
% Continuity Equation
dU1_dt_predictor(j) = -dF1_dx;
% Momentum Equation
dU2_dt_predictor(j) = -dF2_dx + J2_initial;
% Energy Equation
dU3_dt_predictor(j) = -dF3_dx;
% Performing the solution update
U1_initial(j) = U1_initial(j) + dU1_dt_predictor(j)*dt;
U2_initial(j) = U2_initial(j) + dU2_dt_predictor(j)*dt;
U3_initial(j) = U3_initial(j) + dU3_dt_predictor(j)*dt;
end
%
% Calculating the new F1 , F2 and F3 to be used in the corrector method
%
F1_initial = U2_initial;
F2_initial = ( U2_initial.^2 )./U1_initial + (gamma-1)/gamma .* ...
( U3_initial- (gamma/2) * (U2_initial.^2)./U1_initial );
F3_initial = U2_initial .* U3_initial ./ U1_initial + ...
(gamma-1) .* U2_initial.* ( U3_initial ./ U1_initial - ...
gamma/2 * U2_initial.^2 ./ U1_initial.^2 );
%
% Now we have calculated all the paramters after performing the
% predictor step
%
% Now after the predictor step we shall be performing the Corrector
% step
%
% Corrector step uses the backward differencing scheme for its spatial
% terms
%
% In the calculation of the spatial terms we will be using the values
% we calculated from the predictor step
for j = 2:n-1
% Calculating the individual terms in the equations so as to easy
% to program and easy to debug
%
% The spactial terms using the Backward differencing
dF1_dx = ( F1_initial(j) - F1_initial(j-1) )/dx;
dF2_dx = ( F2_initial(j) - F2_initial(j-1) )/dx;
dF3_dx = ( F3_initial(j) - F3_initial(j-1) )/dx;
dlogA_dx = ( log(A(j)) - log(A(j-1)) )/dx;
J2_initial = (gamma-1)/(gamma)*( U3_initial(j) - (gamma/2) *...
U2_initial(j)^2/U1_initial(j) ) * dlogA_dx ;
% Performing the predictor method
% Continuity Equation
dU1_dt_corrector(j) = -dF1_dx;
% Momentum Equation
dU2_dt_corrector(j) = -dF2_dx + J2_initial;
% Energy Equation
dU3_dt_corrector(j) = -dF3_dx;
end
%
% Now after performing the corrector step we will be finding the
% average derivative
dU1_dt_average = 0.5 * ( dU1_dt_corrector + dU1_dt_predictor );
dU2_dt_average = 0.5 * ( dU2_dt_corrector + dU2_dt_predictor );
dU3_dt_average = 0.5 * ( dU3_dt_corrector + dU3_dt_predictor );
% Performing the final solution update
for j = 2:n-1
% Updating the values of the initial variable
U1_initial(j) = U1_old(j) + dU1_dt_average(j)*dt;
U2_initial(j) = U2_old(j) + dU2_dt_average(j)*dt;
U3_initial(j) = U3_old(j) + dU3_dt_average(j)*dt;
end
% Applying the boundary conditions
%
% At inlet
U2_initial(1) = 2*U2_initial(2) - U2_initial(3);
v_first_element = U2_initial(1)./U1_initial(1);
U3_initial(1) = U1_initial(1)*( 1/(gamma-1) +...
(gamma/2)*v_first_element^2 );
% At the outlet boundary
U1_initial(n) = 2*U1_initial(n-1) - U1_initial(n-2);
U2_initial(n) = 2*U2_initial(n-1) - U2_initial(n-2);
U3_initial(n) = 2*U3_initial(n-1) - U3_initial(n-2);
% Finding the final variables
rho_initial = U1_initial./A;
v_initial = U2_initial./U1_initial;
T_initial = (gamma-1).*( U3_initial./U1_initial - ...
(gamma/2)*v_initial.^2 );
P_initial = rho_initial.*T_initial;
% Calculating the flux terms
F1_initial = U2_initial;
F2_initial = ( U2_initial.^2 )./U1_initial + (gamma-1)/gamma .* ...
( U3_initial- (gamma/2) * (U2_initial.^2)./U1_initial );
F3_initial = U2_initial .* U3_initial ./ U1_initial + ...
(gamma-1) .* U2_initial.* ( U3_initial ./ U1_initial - ...
gamma/2 * U2_initial.^2 ./ U1_initial.^2 );
% Finding the state error to ensure convergence into steady state
error_U1 = max(max(abs( U1_initial - U1_old )));
error_U2 = max(max(abs( U2_initial - U2_old )));
error_U3 = max(max(abs( U3_initial - U3_old )));
% Finding the maximum error amongst them
state_error = max([ error_U1 ; error_U2 ; error_U3 ]);
% Calculating the normalised mass-flow rate
mass_rate_normalized_conservative = rho_initial.*A.*v_initial;
% Finding the normalised variables at the nozzle cross section
r(nt) = rho_initial(point);
m(nt) = mass_rate_normalized_conservative(point);
T(nt) = T_initial(point);
P(nt) = P_initial(point);
v(nt) = v_initial(point);
% Plotting results at an interval of 20 time steps
% Writing the frames only in the multiple of 10
if rem(nt,a) == 0
% Updating the loop counter to write movies
mm = mm + 1;
% Plotting some results as a continuous animation
figure(2);clf;
% Resizing the window
set(gcf, 'Position', [100, 100, 1400, 1200]);
% subplot 1 - Velocity vs. Length
subplot(221)
plot(xx,v_initial,'linew',2,'color','b');
xlabel('Non-Dimensional X-Axis');
ylabel({'Non-Dimensional';' Veloocity'});
title(['Velocity vs. Length at Time = ' num2str(nt*dt)...
' Seconds']);
set(gca,'fontsize',17);
grid on;
grid minor;
% subplot 2 - Density vs. Length
subplot(222)
plot(xx,rho_initial,'linew',2,'color','r');
xlabel('Non-Dimensional X-Axis');
ylabel({'Non-Dimensional';' Density'});
title(['Density vs. Length at Time = ' num2str(nt*dt)...
' Seconds']);
set(gca,'fontsize',17);
grid on;
grid minor;
% subplot 3 - Temperature vs. Length
subplot(223)
plot(xx,T_initial,'linew',2,'color','m');
xlabel('Non-Dimensional X-Axis');
ylabel({'Non-Dimensional';' Temperature'});
title(['Temperature vs. Length at Time = ' num2str(nt*dt) ...
' Seconds']);
set(gca,'fontsize',17);
grid on;
grid minor;
% subplot 4 - Pressure vs. Length
subplot(224)
plot(xx,P_initial,'linew',2,'color','g');
xlabel('Non-Dimensional X-Axis');
ylabel({'Non-Dimensional';' Pressure'});
title(['Pressure vs. Length at Time = ' num2str(nt*dt) ...
' Seconds']);
set(gca,'fontsize',17);
grid on;
grid minor;
% update movie frame
M_1(mm) = getframe(gcf);
% Now mass-flow rate normalised
figure(3);clf;
set(gcf,'Position',[ 0 0 1400 1200 ]);
plot(xx,mass_rate_normalized_conservative,'linew',2,'color','r');
xlabel('Non-Dimensional X-Axis');
ylabel({'Non-Dimensional';' Mass-Flow Rate'});
ylim([0 1]);
title(['Mass-Flow Rate vs. Length at Time = ' num2str(nt*dt) ...
' Seconds']);
set(gca,'fontsize',17);
grid on;
grid minor;
M_2(mm) = getframe(gcf);
% Pausing to give animation effect
pause(0.02);
end
end
% Assigning the final variables
iterations_to_steady = nt;
% The primitive variables
m_rate = mass_rate_normalized_conservative;
% The primitive variables
rho_final = rho_initial;
v_final = v_initial;
T_final = T_initial;
P_final = P_initial;
% Time taken to solve the case
Time_taken = nt*dt;
endThe function is made highly modular leaving up to the user to decide the rate of frames for animation, provide the tolerance required for the steady-state conditions and returning the Final Time taken for the simulation and returning the variation of the primitive variables at a specified grid point with a number of iterations. With the help of following code, the steady-state information is obtained.
%% Computing the iterations using the solver
% I will implement the solver I built to solve the conservative form of the
% MAcCormack Method
[rho_final,v_final,T_final,P_final,m_rate,r,v,T,P,m,...
iterations_to_steady,M_1,M_2,Time_taken] = ...
conservative_MacCormack_solver(gamma,rho_initial,...
T_initial,A,P_initial,tol,xx,n,interval_movie_plot,point);
% Assigining the number of Iterations taken to nt
nt = iterations_to_steady;
% Printing the message for the user to enter anything
fprintf('Please Press Spacebar to continue : ');
fprintf('\n');
pause;
% Printing the final message for the convergence time
fprintf(['It took ' num2str(Time_taken) ...
' Seconds for the solution to reach steady state !! ']);
fprintf('\n');
fprintf('Please Press Spacebar to continue : ');
fprintf('\n');
pause;
The following output is obtained:
And now after the plotting it is time to make a movie and plot the final plots at steady state which is done using the following command:
%% Final plots
% Plotting some final results
figure(4);clf;
% Resizing the window
set(gcf, 'Position', [100, 100, 1400, 1200]);
% subplot 1 - Velocity vs. Length
subplot(221)
plot(xx,v_final,'linew',2,'color','b');
xlabel('Non-Dimensional X-Axis');
ylabel({'Non-Dimensional';' Veloocity'});
title(['Velocity vs. Length at Time = ' num2str(Time_taken)...
' Seconds']);
set(gca,'fontsize',17);
grid on;
% subplot 2 - Density vs. Length
subplot(222)
plot(xx,rho_final,'linew',2,'color','r');
xlabel('Non-Dimensional X-Axis');
ylabel({'Non-Dimensional';' Density'});
title(['Density vs. Length at Time = ' num2str(Time_taken)...
' Seconds']);
set(gca,'fontsize',17);
grid on;
% subplot 3 - Temperature vs. Length
subplot(223)
plot(xx,T_final,'linew',2,'color','m');
xlabel('Non-Dimensional X-Axis');
ylabel({'Non-Dimensional';' Temperature'});
title(['Temperature vs. Length at Time = ' num2str(Time_taken) ...
' Seconds']);
set(gca,'fontsize',17);
grid on;
% subplot 4 - Pressure vs. Length
subplot(224)
plot(xx,P_final,'linew',2,'color','g');
xlabel('Non-Dimensional X-Axis');
ylabel({'Non-Dimensional';' Pressure'});
title(['Pressure vs. Length at Time = ' num2str(Time_taken) ...
' Seconds']);
set(gca,'fontsize',17);
grid on;
% Now mass-flow rate normalised
figure(5);clf;
set(gcf,'Position',[ 0 0 1400 1200 ]);
plot(xx,m_rate,'linew',2,'color','r');
xlabel('Non-Dimensional X-Axis');
ylabel({'Non-Dimensional';' Mass-Flow Rate'});
ylim([0 1]);
title(['Mormalized Mass-Flow Rate vs. Length at Time = '...
num2str(Time_taken) ' Seconds']);
set(gca,'fontsize',17);
grid on;
grid minor;The following plots are obtained:
As it can be seen the mass rate is quite smooth globally and is near to 0.59 which agrees with the theoretical calculations. But if we see zoomed in mass-flow rate we see local instabilities, this occurs mainly due to the discretization error and the rounding off error. The following is the plot for mass-flow rate but a bit zoomed in:
As can be seen there are local instabilities.
To obtain the animation of the simulation the following code is implemented :
%% Making an Animation
% Writing out a movie
% For the variatio of parameters
videofile1 = VideoWriter('Primitive Parameters conservative MacCormack.avi');
videofile1.FrameRate = 5;
open(videofile1)
writeVideo(videofile1,M_1);
close(videofile1);
% For the varition of mass-flow rate
videofile2 = VideoWriter('Mass-Flow normalized rate conservative MacCormack.avi');
videofile2.FrameRate = 5;
open(videofile2)
writeVideo(videofile2,M_2);
close(videofile2);The final animation was obtained for the animation of the state variables:
The final animation was obtained for the animation of the normalized mass-flow rate:
It is clearly evident that at the initial time steps the fluctuations in the mass-rate overall is very high. This is because of the fact of the presence of the derivative terms and their counterparts which take time to damp. As can be seen from the animation the solution becomes nearly stable at 10 seconds. The extra simulation is to ensure that the solution has indeed converged to a specified tolerance by the user.
Some more plots at the intermediate time steps:
It is clearly seen that the variations in the variables after 10 seconds almost seem to be negligible. Here we can say that the solution has reached a well enough steady-state for the human eye, however for machine precision and accuracy it continues to 41 seconds in order the allow the consecutive iterates not to differ by more than the specified tolerance.
We must also look at the plots of the variation of the primitive variables at a specific grid point per iteration which is implemented with the help of the following MATLAB Code after running the upper codes.
%% Plotting the variations of non-dimensional variables at grid point i
figure(6);clf;
set(gcf,'Position',[ 0 0 1400 1200 ]);
plot(m,'linew',2,'color','r');
xlabel('Iterations');
ylabel({'Non-Dimensional';' Mass-Flow Rate'});
title(['Normalized MASS-RATE vs. Iterations at grid point ' ...
num2str(point)]);
set(gca,'fontsize',17);
grid on;
grid minor;
figure(7);clf;
set(gcf,'Position',[ 0 0 1400 1200 ]);
plot(P,'linew',2,'color','r');
xlabel('Iterations');
ylabel({'Non-Dimensional';' Temperature'});
title(['Normalized Temperature vs. Iterations at grid point ' ...
num2str(point)]);
set(gca,'fontsize',17);
grid on;
grid minor;
figure(8);clf;
set(gcf,'Position',[ 0 0 1400 1200 ]);
plot(r,'linew',2,'color','r');
xlabel('Iterations');
ylabel({'Non-Dimensional';' Density'});
title(['Normalized Density vs. Iterations at grid point ' ...
num2str(point)]);
set(gca,'fontsize',17);
grid on;
grid minor;
figure(9);clf;
set(gcf,'Position',[ 0 0 1400 1200 ]);
plot(P,'linew',2,'color','r');
xlabel('Iterations');
ylabel({'Non-Dimensional';' Pressure'});
title(['Normalized Pressure vs. Iterations at grid point ' ...
num2str(point)]);
set(gca,'fontsize',17);
grid on;
grid minor;
The following graphs were obtained :
Important Observations:
1.) Early times, the time derivatives are large and hence they oscillate a lot. But it is seen from the above graphs that after some iterations they settle down and are not changing much less which is exactly what we want for our numerical scheme i.e. to converge to a solution.
2.) We can see that our numerical scheme indeed produces accurate enough results as the mass rate initially assumed to be 0.59 is achieved at the steady-state which is visually confirmed from the last graph when we see the curve settling down to a value of 0.59.
The above simulation was performed with the help of 31 grid points. Now we shall change the points to 51 and 61 and run the iteration for 1400 time steps and shall see the results:
1.) For n = 51 points
The following results were obtained by running it for 1400 time loop iterations:
It is seen that for n= 51 grid points we get very good accuracy in predicting the results and which makes sense also as higher the grid points the more accuracy one tends to get in the numerical scheme.
2.) For n = 61 points
The following animation was obtained for the variation of the mass-flow rate with time
link : https://youtu.be/Dip_dd5wXxc
And for the primitive variables :
link - https://youtu.be/2J3UzW1glww
As it is evident from the animation there are instabilities at the end of the boundaries which causes the whole solution to blow up. This actually has to do with the selection of the number of grid points. It may seem that as we increase the number of grid points we get more accuracy. However, it is not the only factor in determining the accuracy and stability of a numerical scheme. There is something called a CFL constant which governs the stability of the solution. See the plots that were obtained when the solution started to blow up.
The CFL condition expresses that the distance that any information travels during the timestep length within the mesh must be lower than the distance between mesh elements. In other words, information from a given cell or mesh element must propagate only to its immediate neighbors.
It is evident from the data and the simulations that the CFL constant plays a very crucial role in determining the stability of the solution. This CFL constant is directly dependent upon the grid size `Delta x` and the time step size `Delta t`. The selection of these parameters should be done with care if we allow the convergence of the solution otherwise the solution will blow up. In general, the notion is that the CFL constant is bound by a constant which in case of the explicit finite difference methods is 1. If the CFL constant is greater than 1 then the solution is unstable and will blow up. If the CFL number is well within the bounds of this constant then the solution is stable and the scheme will converge. If the grid sizes are made to account for more accuracy proportional changes must be made in the time step size to account for the convergence and the stability.
It has to with the fact that how fast the information travels across the grid point. As the grid points are increased the effective grid spacing is decreased and hence a proportionate increase in the time step size should be made in order to compensate for the effect to allow for the convergence of the result.
The CFL condition is commonly confused with linear stability conditions or nonlinear stability conditions. However, it’s important to emphasize that it is not a sufficient condition for stability, and other stability conditions are generally more restrictive than the CFL condition. It is therefore imperative to choose the CFL constant wisely in order for the simulation to give meaningful results.
Non-Conservative Approach
Now with a basic knowledge of the governing equations and the MacCormack's Technique we shall move our study to solving the 1-D nozzle flow problem using the non-conservative approach. Nut first it is important to derive the governing equations in the non-conservative form. These are derived by segregating each and every term in the conservative equations. The following were the conservative equations derived :
`frac(partial)(partialt)(rhoA) + frac(partial)(partialx)(rhoAv) = 0`
`frac(partial)(partialt)(rhovA) + frac(partial)(partialx)(rhov^2A) = -Afrac(partialP)(partialx)`
`frac(partial)(partialt)(rho[e + frac(v^2)(2)]A) + frac(partial)(partialx)(rhoAv[e + frac(v^2)(2)])= -frac(partial)(partialx)(PAv)`
The Continuity equation :
`frac(partial)(partialt)(rhoA) + frac(partial)(partialx)(rhoAv) = 0`
As the Area vector is independent of time ( assuming no deformations in the shape of the nozzle ) we can write the equation as :
`Afrac(partial)(partialt)(rho) +Av frac(partial)(partialx)(rho) + vrhofrac(partial)(partialx)(A) + Arhofrac(partial)(partialx)(v) = 0`
We can simplify it as :
`frac(partial)(partialt)(rho) +v frac(partial)(partialx)(rho) + 1/Avrhofrac(partial)(partialx)(A) + rhofrac(partial)(partialx)(v) = 0`
Taking all the spatial terms to the right hand side we can write this as :
`frac(partial)(partialt)(rho)= -v frac(partial)(partialx)(rho) - 1/Avrhofrac(partial)(partialx)(A) - rhofrac(partial)(partialx)(v)`
which is the non-conservative form of the continuity equation.
The Momentum Equation
`frac(partial)(partialt)(rhovA) + frac(partial)(partialx)(rhov^2A) = -Afrac(partialP)(partialx)`
We can expand the time and spatial derivative terms as:
`vfrac(partial)(partialt)(rhoA)+Arhofrac(partial)(partialt)(v) + vfrac(partial)(partialx)(rhovA) +rhovAfrac(partial)(partialx)(v) = -Afrac(partialP)(partialx)`
If one notices carefully :
`vfrac(partial)(partialt)(rhoA)+ vfrac(partial)(partialx)(rhovA) = v [ frac(partial)(partialt)(rhoA)+ frac(partial)(partialx)(rhovA)] = v xx 0 = 0` from the continuity equation, hence the equation simplifies to. :
`Arhofrac(partial)(partialt)(v) +rhovAfrac(partial)(partialx)(v) = -Afrac(partialP)(partialx)`
or we can segregate the time and the spatial derivative terms as :
`Arhofrac(partial)(partialt)(v) =-rhovAfrac(partial)(partialx)(v) -Afrac(partialP)(partialx)`
Dividing throughout by A we get :
`rhofrac(partial)(partialt)(v) =-rhovfrac(partial)(partialx)(v) -frac(partialP)(partialx)`
which is the non-conservative form of the momentum Equation.
The Energy Equation
`frac(partial)(partialt)(rho[e + frac(v^2)(2)]A) + frac(partial)(partialx)(rhoAv[e + frac(v^2)(2)])= -frac(partial)(partialx)(PAv)`
This equation can be simplified with the help of momentum equation as:
`frac(partial)(partialt)(rhoeA) + frac(partial)(partialx)(rhoevA) = -Pfrac(partial)(partialx)(Av)`
`rhoAfrac(partial)(partialt)(e) + eAfrac(partial)(partialt)(rho) + efrac(partial)(partialx)(rhoAv) + rhoAvfrac(partial)(partialx)(e) = -P[vfrac(partial)(partialx)(A) + Afrac(partial)(partialx)(v)]`
It can be noted carefully that `Aefrac(partial)(partialt)(rho) + efrac(partial)(partialt)(rhoAv) = e [ frac(partial)(partialt)(Arho)+frac(partial)(partialt)(rhoAv)]=exx0 = 0`. So we can write the equation as:
`rhoAfrac(partiale)(partialt) + rhoAvfrac(partiale)(partialx)= - rhoRT[vfrac(partial)(partialx)(A) + Afrac(partial)(partialx)(v)]`
Dividing throughout by A we get
`rhofrac(partiale)(partialt) + rhovfrac(partiale)(partialx) =- rhoRT[v/Afrac(partial)(partialx)(A) + frac(partial)(partialx)(v)]`
Using the substitution `e = C_vT` we can write the equation as:
`(partialT)/(partialt) + v(partialT)/(partialx) = -(gamma-1)T[(partialv)/(partialx) + v (partial(lnA))/(partialx)]`
Segregating the time and spatial derivative terms we can write:
`(partialT)/(partialt) =- v(partialT)/(partialx) -(gamma-1)T[(partialv)/(partialx) + v (partial(lnA))/(partialx)]`
which is the non-conservative form of the Energy Equation.
Now we shall non-dimensionalize the equations for their easy implementation.
We choose the non-dimensional parameters as :
- `x' = frac(x)(L)` where L is the length of the nozzle
- `t' = frac(t)(L//v_0)` where `v_0` is the sonic speed
- `rho' = rho/rho_0`
- `T' = T/T_0`
- `v' = v/v_0` and finally
- `A' = A/A_0`
where the values `rho_0, T_0,v_0 and A_0` are the values at the throat of the nozzle.
Substituting these variables in the continuity equation we get:
`frac(partial)(partialt')(rho')= -v' frac(partial)(partialx')(rho') - v'rho'frac(partial)(partialx')(lnA') - rho'frac(partial)(partialx')(v')`
Substituting these variables in the momentum equation we get:
`rho'frac(partial)(partialt')(v') =-rho'v'frac(partial)(partialx')(v') -1/gammafrac(partialP')(partialx')`
Using the fact that `P' = rho'T'`
`rho'frac(partial)(partialt')(v') =-rho'v'frac(partial)(partialx')(v') -1/gamma[rho'frac(partialT')(partialx') + T'frac(partialrho')(partialx')]`
which is the non-dimensional form of the non-conservative momentum equation.
Substituting these variables in the energy equation we get:
`(partialT')/(partialt') =- v(partialT')/(partialx') -(gamma-1)T'[(partialv')/(partialx') + v '(partial(lnA'))/(partialx')]`
which is the non-dimensional form of the non-conservative energy equation.
So now we have overall three non-conservative equations as follows:
`frac(partial)(partialt')(rho')= -v' frac(partial)(partialx')(rho') - v'rho'frac(partial)(partialx')(lnA') - rho'frac(partial)(partialx')(v')`
`rho'frac(partial)(partialt')(v') =-rho'v'frac(partial)(partialx')(v') -1/gamma[rho'frac(partialT')(partialx') + T'frac(partialrho')(partialx')]`
and
`(partialT')/(partialt') =- v(partialT')/(partialx') -(gamma-1)T'[(partialv')/(partialx') + v '(partial(lnA'))/(partialx')]`
A very important observation to note here is the fact that the time derivative terms involve primitive variables individually and hence we need not make any substitution into the equations to modify them.
Initial Conditions Setup
The Initial conditions are chosen as:
`A = 1 + 2.2(x-1.5)^2`
`rho = 1 - 0.3146x`
`T = 1 - 0.2314x` and
`v = (0.1 + 1.09x)T^0.5`
The Following MATLAB Code is implemented to initiate the initial conditions and the number of grid points, movie interval plot, etc just as in conservative approach:
clear;
close all;
clc;
%
% In this program I shall be coding the nozzle problem in 1-D using the
% MacCormack's method with the use of the predictor corrector method
%
% The problem domain is selected as [0,3] for the non-dimensional variable
% x
%
% The number of grid points used for the discretization including the
% boundary points are n = 31
%
% The initial temperature profiles are assumed @t=0
%
% The equations gorverning the nozzle flow are non-linear PDE and are
% non-conservative
%
% Defining the number of grid points
n = 31;
% Defining the problem domain
xx = linspace(0,3,n);
% Defining the grid size
dx = xx(2) - xx(1);
% Define the gamma
gamma = 1.4;
% Defining the Area profile
A = 1 + 2.2*(xx-1.5).^2;
% Defining the initial conditions @(t=0)
rho_initial = 1 - 0.3146*xx;
T_initial = 1 - 0.2314*xx;
v_initial = (0.1 + 1.09*xx).*T_initial.^0.5;
% We can define the initial pressure
%
% P = rho * T in the non-dimensional form
P_initial = rho_initial.*T_initial;
% Defining the tolerance
tol = 1e-6;
% Defining the time interval for movie plot
interval_movie_plot = 10;
% Defining the point at which we want to study the variation
point = round(n/6);
The following Initial plots are obtained :
%% Plotting the initial conditions
figure(1);clf;
set(gcf,'Position',[0 0 1600 1600]);
subplot(221);
plot(xx,rho_initial,'r-','linew',2);
xlabel('Non-Dimensional X-Axis');
ylabel({'Non-Dimensional';' Density'});
title('Initial Density');
set(gca,'fontsize',17);
grid on;
grid minor;
subplot(222);
plot(xx,v_initial,'r-','linew',2);
xlabel('Non-Dimensional X-Axis');
ylabel({'Non-Dimensional';' Velocity'});
title('Initial Velocity');
set(gca,'fontsize',17);
grid on;
grid minor;
subplot(223);
plot(xx,T_initial,'r-','linew',2);
xlabel('Non-Dimensional X-Axis');
ylabel({'Non-Dimensional';' Temperature'});
title('Initial Temperature');
set(gca,'fontsize',17);
grid on;
grid minor;
subplot(224);
plot(xx,P_initial,'r-','linew',2);
xlabel('Non-Dimensional X-Axis');
ylabel({'Non-Dimensional';' Pressure'});
title('Initial Density');
set(gca,'fontsize',17);
grid on;
grid minor;A solver is built in MATLAB which takes the initial condition from the above code and solves the quasi 1-D nozzle flow problem using the non-conservative approach. The code implemented to construct the solver is depicted below:
function [rho_final,v_final,T_final,P_final,m_rate,r,v,T,P,m,...
iterations_to_steady,M_1,M_2,Time_taken] = ...
non_conservative_MacCormack_solver(gamma,...
rho_initial,T_initial,A,P_initial,v_initial,tol,xx,n,a,point)
%
% This function will take input as the initial density and initial
% Temperature as the input and the steady state tolerance desired and will
% produce M_1 and M_2 movie loop counters holding the frames, the primitive
% state variables at the final step and the variation of non-dimensional
% variables with iterations at the grid point specified by the user
%
% The non-conservative equations of the 1-D quasi flow are given by:
%
% d/dt(rho) = .....
% d/dt(v) = .....
% d/dt(T) = ......
%
% Basically in the non-conservative form of the governing equation the
% primitive variables appear individually and hence we don't need to make
% any fancy sustitution to solve them
%
% Defining the grid spacing
dx = xx(2) - xx(1);
% Defining the time step
dt = 0.01;
% Initialising the loop counter
nt = 0;
% Defining the movie loop counter
mm = 0;
% Defining the state error initially
state_error = 9e9;
%
% It is imperative to note that these variables are non-dimensional
%
% Now we begin our MacCormack method
%
% Main outer time loop
while(state_error>tol)
% Increement the loop counter
nt = nt + 1;
% Keep a copy of the old variables
rho_old = rho_initial;
v_old = v_initial;
T_old = T_initial;
P_old = P_initial;
% Predictor Method
%
% First we shall calculate the predictor variables using the forward
% differencing scheme
for j = 2:n-1
% Calculating the individual terms in the equations so as to easy
% to program and easy to debug
%
% The spactial terms using the forward differencing
drho_dx = ( rho_initial(j+1) - rho_initial(j) )/dx;
dv_dx = ( v_initial(j+1) - v_initial(j) )/dx;
dlogA_dx = ( log(A(j+1)) - log(A(j)) )/dx;
dT_dx = ( T_initial(j+1) - T_initial(j) )/dx;
% The main time derivative terms
%
% Continuity Equation
drho_dt_predictor(j) = -v_initial(j)*drho_dx - rho_initial(j)*dv_dx - ...
rho_initial(j)*v_initial(j)*dlogA_dx;
% Momentum Equation
dv_dt_predictor(j) = -v_initial(j)*dv_dx - 1/(gamma) * ...
( dT_dx + (T_initial(j)./rho_initial(j))* drho_dx );
% Energy Equation
dT_dt_predictor(j) = -v_initial(j)*dT_dx -...
(gamma-1) * T_initial(j) * ( dv_dx + v_initial(j)*dlogA_dx );
% Updating the variables in the predictor step:
rho_initial(j) = rho_initial(j) + drho_dt_predictor(j)*dt;
v_initial(j) = v_initial(j) + dv_dt_predictor(j)*dt;
T_initial(j) = T_initial(j) + dT_dt_predictor(j)*dt;
end
% Now after the predictor step we shall be performing the Corrector
% step
%
% Corrector step uses the backward differencing scheme for its spatial
% terms
%
% In the calculation of the spatial terms we will be using the values
% we calculated from the predictor step
for j = 2:n-1
% Calculating the individual terms in the equations so as to easy
% to program and easy to debug
%
% The spactial terms using the forward differencing
drho_dx = ( rho_initial(j) - rho_initial(j-1) )/dx;
dv_dx = ( v_initial(j) - v_initial(j-1) )/dx;
dlogA_dx = ( log(A(j)) - log(A(j-1)) )/dx;
dT_dx = ( T_initial(j) - T_initial(j-1) )/dx;
% The main time derivative terms
%
% Continuity Equation
drho_dt_corrector(j) = -v_initial(j)*drho_dx - rho_initial(j)*dv_dx - ...
rho_initial(j)*v_initial(j)*dlogA_dx;
% Momentum Equation
dv_dt_corrector(j) = -v_initial(j)*dv_dx - 1/(gamma) * ...
( dT_dx + (T_initial(j)/rho_initial(j))* drho_dx );
% Energy Equation
dT_dt_corrector(j) = -v_initial(j)*dT_dx - (gamma-1)*T_initial(j)*...
( dv_dx + v_initial(j)*dlogA_dx );
end
% Now we have performed the corrector step successfully
%
% Now we will calculate the average derivatives
drho_dt_average = 0.5*(drho_dt_corrector + drho_dt_predictor);
dv_dt_average = 0.5*(dv_dt_corrector + dv_dt_predictor);
dT_dt_average = 0.5*(dT_dt_corrector + dT_dt_predictor );
% Updating the final values at the end of 1 time loop iteration
for j = 2:n-1
% Updating the final variables
v_initial(j) = v_old(j) + dv_dt_average(j)*dt;
rho_initial(j) = rho_old(j) + drho_dt_average(j)*dt;
T_initial(j) = T_old(j) + dT_dt_average(j)*dt;
end
% Updating the boundary conditions
v_initial(1) = 2*v_initial(2) - v_initial(3);
% The boundary conditions for T and rho are the same at the beginning
% point
% End boundary conditions
v_initial(n) = 2*v_initial(n-1) - v_initial(n-2);
rho_initial(n) = 2*rho_initial(n-1) - rho_initial(n-2);
T_initial(n) = 2*T_initial(n-1) - T_initial(n-2);
P_initial = rho_initial.*T_initial;
% Calculating the normalised mass-flow rate
mass_rate_normalized_non_conservative = rho_initial.*A.*v_initial;
% Defining the differences between the previous and the current time
% iteration values
error_v = max(max(abs(v_initial-v_old)));
error_rho = max(max(abs(rho_initial-rho_old)));
error_T = max(max(abs(T_initial-T_old)));
error_P = max(max(abs(P_initial-P_old)));
%
% Maximum of these will correspond to the state error
state_error = max( [ error_v ; error_rho ; error_T ; error_P ] );
% Returning the primitive state variables with variation of iterations
r(nt) = rho_initial(point);
v(nt) = v_initial(point);
T(nt) = T_initial(point);
P(nt) = P_initial(point);
m(nt) = mass_rate_normalized_non_conservative(point);
% Now
% Writing the frames only in the multiple of 10
if rem(nt,a) == 0
% Updating the loop counter to write movies
mm = mm + 1;
% Plotting some results as a continuous animation
figure(2);clf;
% Resizing the window
set(gcf, 'Position', [100, 100, 1400, 1200]);
% subplot 1 - Velocity vs. Length
subplot(221)
plot(xx,v_initial,'linew',2,'color','b');
xlabel('Non-Dimensional X-Axis');
ylabel({'Non-Dimensional';' Velocity'});
title(['Velocity vs. Length at Time = ' num2str(nt*dt)...
' Seconds']);
set(gca,'fontsize',17);
grid on;
% subplot 2 - Density vs. Length
subplot(222)
plot(xx,rho_initial,'linew',2,'color','r');
xlabel('Non-Dimensional X-Axis');
ylabel({'Non-Dimensional';' Density'});
title(['Density vs. Length at Time = ' num2str(nt*dt)...
' Seconds']);
set(gca,'fontsize',17);
grid on;
% subplot 3 - Temperature vs. Length
subplot(223)
plot(xx,T_initial,'linew',2,'color','m');
xlabel('Non-Dimensional X-Axis');
ylabel({'Non-Dimensional';' Temperature'});
title(['Temperature vs. Length at Time = ' num2str(nt*dt) ...
' Seconds']);
set(gca,'fontsize',17);
grid on;
% subplot 4 - Pressure vs. Length
subplot(224)
plot(xx,P_initial,'linew',2,'color','g');
xlabel('Non-Dimensional X-Axis');
ylabel({'Non-Dimensional';' Pressure'});
title(['Pressure vs. Length at Time = ' num2str(nt*dt) ...
' Seconds']);
set(gca,'fontsize',17);
grid on;
% update movie frame
M_1(mm) = getframe(gcf);
% Now mass-flow rate normalised
figure(3);clf;
set(gcf,'Position',[ 0 0 1400 1200 ]);
plot(xx,mass_rate_normalized_non_conservative,'linew',2,'color','r');
xlabel('Non-Dimensional X-Axis');
ylabel({'Non-Dimensional';' Mass-Flow Rate'});
title(['Mass-Flow Rate vs. Length at Time = ' num2str(nt*dt) ...
' Seconds']);
set(gca,'fontsize',17);
grid on;
M_2(mm) = getframe(gcf);
% Pausing to give animation effect
pause(0.01);
end
end
% Returning the final variables
m_rate = mass_rate_normalized_non_conservative;
rho_final = rho_initial;
v_final = v_initial;
T_final = T_initial;
P_final = P_initial;
% Returning the number of iterations taken to reach steady state
iterations_to_steady = nt;
% Returning the time taken for simulation
Time_taken = nt*dt;
end
The function is made highly modular leaving upto the user to decide the rate of frames for animation, provide the tolerance required for the steady-state conditions and returning the Final Time taken for the simulation and returning the variation of the primitive variables at a specified grid point with number of iterations. With the help of following code the steady state information is obtained.
%% Computing the iterations using the solver
% I will implement the solver I built to solve the conservative form of the
% MAcCormack Method
[rho_final,v_final,T_final,P_final,m_rate,r,v,T,P,m,...
iterations_to_steady,M_1,M_2,Time_taken] = ...
non_conservative_MacCormack_solver(gamma,...
rho_initial,T_initial,A,P_initial,v_initial...
,tol,xx,n,interval_movie_plot,point);
% Pausing for the user to press any key
fprintf('\n');
fprintf('Pause any key to Continue for the final plots');
fprintf('\n');
pause;
% Printing the message out to the screen
fprintf(['It took ' num2str(iterations_to_steady) ' Iterations in the time Loop' ...
' or ' num2str(Time_taken) ' Seconds for the simulation to converge']);
fprintf('\n');
fprintf('Pause any key to Continue for the final plots');
fprintf('\n');
pause;
The following output is obtained:
And now after the plotting, it is time to make a movie and plot the final plots at steady state which is done using the following command:
%% Final Plots
%
% Plotting final results
figure(4);clf;
% Resizing the window
set(gcf, 'Position', [100, 100, 1400, 1200]);
% subplot 1 - Velocity vs. Length
subplot(221)
plot(xx,v_final,'linew',2,'color','b');
xlabel('Non-Dimensional X-Axis');
ylabel({'Non-Dimensional';' Velocity'});
title(['Velocity vs. Length at Time = ' num2str(Time_taken)...
' Seconds']);
set(gca,'fontsize',17);
grid on;
% subplot 2 - Density vs. Length
subplot(222)
plot(xx,rho_final,'linew',2,'color','r');
xlabel('Non-Dimensional X-Axis');
ylabel({'Non-Dimensional';' Density'});
title(['Density vs. Length at Time = ' num2str(Time_taken)...
' Seconds']);
set(gca,'fontsize',17);
grid on;
% subplot 3 - Temperature vs. Length
subplot(223)
plot(xx,T_final,'linew',2,'color','m');
xlabel('Non-Dimensional X-Axis');
ylabel({'Non-Dimensional';' Temperature'});
title(['Temperature vs. Length at Time = ' num2str(Time_taken) ...
' Seconds']);
set(gca,'fontsize',17);
grid on;
% subplot 4 - Pressure vs. Length
subplot(224)
plot(xx,P_final,'linew',2,'color','g');
xlabel('Non-Dimensional X-Axis');
ylabel({'Non-Dimensional';' Pressure'});
title(['Pressure vs. Length at Time = ' num2str(Time_taken) ...
' Seconds']);
set(gca,'fontsize',17);
grid on;
figure(5);clf;
% Now mass-flow rate normalised
set(gcf,'Position',[ 0 0 1600 1600 ]);
plot(xx,m_rate,'linew',2,'color','r');
xlabel('Non-Dimensional X-Axis');
ylabel({'Non-Dimensional';' Mass-Flow Rate'});
title(['Mass-Flow Rate vs. Length at Time = ' num2str(Time_taken) ...
' Seconds']);
set(gca,'fontsize',17);
grid on;
The following plots are obtained:
As it can be seen the mass rate is quite smooth globally and is near to 0.59 which agrees with the theoretical calculations. But if we see zoomed in mass-flow rate we see local instabilities, this occurs mainly due to the discretization error and the rounding off error. The following is the plot for mass-flow rate but a bit zoomed in:
As can be seen there are local instabilities.
To obtain the animation of the simulation the following code is implemented :
%% Writing a movie
% setup movie
% Now we can write a movie with this
videofile = VideoWriter('Variation of parameters non-conservative.avi');
% Changing the frame count so as to visualize the angles properly
videofile.FrameRate = 6;
open(videofile);
writeVideo(videofile,M_1);
close(videofile);
% Now we can write a movie with this
videofile2 = VideoWriter('mass_flow_rate variation non-conservative.avi');
% Changing the frame count so as to visualize the angles properly
videofile2.FrameRate = 6;
open(videofile2);
writeVideo(videofile2,M_2);
close(videofile2);
The final animation was obtained for the animation of the state variables:
The final animation was obtained for the animation of the normalized mass-flow rate:
It is clearly evident that at the initial time steps the fluctuations in the mass-rate overall is very high. This is because of the fact of the presence of the derivative terms and their counterparts which take time to damp. As can be seen from the animation the solution becomes nearly stable at 15 seconds. The extra simulation is to ensure that the solution has indeed converged to a specified tolerance by the user.
Some more plots at the intermediate time steps:
It is clearly seen that the variations in the variables after 14 seconds almost seem to be negligible. Here we can say that the solution has reached a well enough steady state for human eye, however for machine precision and accuracy it continues to 17.2 seconds in order the allow the consecutive iterates not to differ by more than the specified tolerance.
We must also look at the plots of the variation of the primitive variables at a specific grid point per iteration which is implemented with the help of the following MATLAB Code after running the upper codes.
%% Variation of primitive variables with iterations
figure(6);clf;
set(gcf,'Position',[ 0 0 1400 1200 ]);
plot(m,'linew',2,'color','r');
xlabel('Iterations');
ylabel({'Non-Dimensional';' Mass-Flow Rate'});
title(['Normalized MASS-RATE vs. Iterations at grid point ' ...
num2str(point)]);
set(gca,'fontsize',17);
grid on;
grid minor;
figure(7);clf;
set(gcf,'Position',[ 0 0 1400 1200 ]);
plot(P,'linew',2,'color','r');
xlabel('Iterations');
ylabel({'Non-Dimensional';' Temperature'});
title(['Normalized Temperature vs. Iterations at grid point ' ...
num2str(point)]);
set(gca,'fontsize',17);
grid on;
grid minor;
figure(8);clf;
set(gcf,'Position',[ 0 0 1400 1200 ]);
plot(r,'linew',2,'color','r');
xlabel('Iterations');
ylabel({'Non-Dimensional';' Density'});
title(['Normalized Density vs. Iterations at grid point ' ...
num2str(point)]);
set(gca,'fontsize',17);
grid on;
grid minor;
figure(9);clf;
set(gcf,'Position',[ 0 0 1400 1200 ]);
plot(P,'linew',2,'color','r');
xlabel('Iterations');
ylabel({'Non-Dimensional';' Pressure'});
title(['Normalized Pressure vs. Iterations at grid point ' ...
num2str(point)]);
set(gca,'fontsize',17);
grid on;
grid minor;The following graphs were obtained :
Important Observations:
1.) Early times, the time derivatives are large and hence they oscillate a lot. But it is seen from the above graphs that after some iterations they settle down and are not changing much less which is exactly what we want for our numerical scheme i.e. to converge to a solution.
2.) We can see that our numerical scheme indeed produces accurate enough results as the mass rate initially assumed to be 0.59 is achieved at the steady-state which is visually confirmed from the last graph when we see the curve settling down to a value of 0.59.
The above simulation was performed with the help of 31 grid points. Now we shall change the points to 51 and 61 and run the iteration for 1400 time steps and shall see the results:
1.) For n = 51 points
The following results were obtained by running it for steady state convergence:
It is seen that for n= 51 grid points we get a very good accuracy in predicting the results and which makes sense also as higher the grid points the more accuracy one tends to get in the numerical scheme.
2.) For n = 131 points
The following graph was obtained for the variation of the mass-flow rate with time:
As is seen from the graph the mass flow rate is quite smooth here as compared to the previous two graphs. It is because of the fact that we have increased the number of grid points and hence we have reduced the numerical discretization error which gives better accuracy in the mass flow rate prediction.
Let us see what happens when we increase the grid points even further.
2.) For n = 101 points
The following animation was obtained for the variation of the mass-flow rate with time
link : https://youtu.be/tt2FPicEjb0
And for the primitive variables :
link - https://youtu.be/4JRNpxcmvuo
As it is evident from the animation there are instabilities at the end of the boundaries which causes the whole solution to blow up. This actually has to do with the selection of number of grid points. It may seem that as we increase the number of grid points we get more accuracy. However, it is not the only factor in determining the accuracy and stability of a numerical scheme. There is something called as a CFL constant which governs the stability of the solution. See the plots that were obtained when the solution started to blow up.
The following graph was obtained for the variation of the mass-flow rate with time:
The following graph was obtained for the variation of the primitive variables with time:
Before concluding the report let us compare the mass-flow rates of the two approaches we just looked at :
Comparison of Mass-Flow rate Using two methods
OBSERVATION:
1.) The conservative form does a better job as to preserving the mass rate as constant in comparison to the non-conservative counterpart which has high fluctuations as compared to the conservative form.
2.) Justification can be given as: in the conservative form, the mass flow rate ( `U_2` ) itself is a dependent parameter in the governing equation. In contrast, the dependent variables in the non-conservative form are the primitive variables, and hence to calculate the mass-flow rate we have to do secondary calculations hence combining the errors from the variables.
CONCLUSION:
- In this report, the study of 1-D quasi nozzle flow is performed. The simulations are carried out using two approaches- conservative and non-conservative forms of the governing equations. Grid independence test is carried out in both the methods to check the effect of the grid size on the stability and the convergence of the solution. The variation of primitive variables with iterations is studied and behavior of the presence of the high derivatives is justified. Lastly, the report is concluded with the comparison of the mass-flow rate that was solved using both the forms, and it is concluded that the conservative approach does a better job in solving the equations.
- However, in the conservative approach, it is a bit mundane to program the code as it involves assuming a lot of variables and recalculating the primitive variables from these variables consuming extra computation power as compared to the non-conservative approach. It is also seen that the time required for the convergence using the conservative approach is considerably higher as compared to the non-conservative approach. One main reason may be attributed to the fact that in a conservative approach the primitive variables are not calculated explicitly and hence it would require much more refinement steps to reach a steady-state solution.
- Also, it is seen that in the conservative approach higher time steps can be taken as compared to the non-conservative counterpart. It is mainly due to the Courant Number ( CFL NUMBER ) which is less for the conservative form allowing it to induce higher time steps.
- Overall, both methods have it own merits and demerits. So hence a proper compromise is necessary in either the accuracy or the simulation time in order to select proper schemes for us to work with.
- The non-conservative part is easy to program and hence should be preferred over the conservative part for the novice programmers with the non-compromisation of reasonable accuracy.