Courses by Software
Courses by Semester
Courses by Domain
Tool-focused Courses
Machine learning
POPULAR COURSES
Post Graduate Program in Hybrid Electric Vehicle Design and AnalysisPost Graduate Program in Computational Fluid DynamicsPost Graduate Program in CADPost Graduate Program in CAEPost Graduate Program in Manufacturing DesignPost Graduate Program in Computational Design and Pre-processingPost Graduate Program in Complete Passenger Car Design & Product DevelopmentSuccess Stories
Analysis and Validation of Quasi One Dimensional De Laval Nozzle using MacCormack Method.
Aim: Simulation of quasi 1D subsonic-supersonic nozzle flow using MacCormack method. Objectives: The objective of this project is to simulate an isentropic flow through a quasi 1D subsonic-supersonic nozzle using MacCormack method for conservative and non-conservative form of governing equations. Theory:…
Tanmay Pathak
updated on 14 Jul 2021
Aim: Simulation of quasi 1D subsonic-supersonic nozzle flow using MacCormack method.
Objectives: The objective of this project is to simulate an isentropic flow through a quasi 1D subsonic-supersonic nozzle using MacCormack method for conservative and non-conservative form of governing equations.
Theory:
The subsonic-supersonic nozzle is also called as "convergent-divergent nozzle" or "de laval nozzle". The steady isentropic flow through a convergent-divergent nozzle is stated below,
Fig (a) : Schematic for subsonic-supersonic isentropic nozzle flow
So here, the flow at the inlet of the nozzle comes from a reservoir where pressure and temperature are denoted as by ρ0ρ0and T0T0 respectively. The cross sectional area of the reservoir is large (theoretically, A→∞A→∞)and hence the velocity is very small (V→0V→0). Thus, ρ0ρ0and T0T0are the stagnation values or total pressure and temperature respectively. The flow expands isentropically to supersonic speeds at the nozzle exits, where the exit pressure, temperature, velocity, and Mach number are denoted by ρeρe,TeTe,VeVeand MeMerespectively. The flow is locally subsonic in the convergent section, sonic at the throat section and supersonic at the divergent section. The sonic flow (M=1M=1) at the throat means that the local velocity at this location is equal to the local speed of sound. Using an asterisk to denote sonic flow values, we have at the throat V=V∗=a∗V=V∗=a∗. Similarly, the sonic flow values of pressure and temperature are denoted by ρ∗ρ∗and T∗T∗, respectively.
Although the area of nozzle changes as a function of distance along the nozzle xx and therefore in reality the flow field is two dimensional (xyxy), we make an assumption that the flow properties vary only with xx, this is tantamount to assume uniform flow properties across any given cross-section. Such flow is defined as Quasi-one-dimensional flow.
In order to solve quasi-one-dimensional problem, we need PDE or Integral form of equations for CFD analysis. These PDE or Integral form of the governing equations are in both conservative and non-conservative form.
Before, writing the governing equations in conservative and non-conservaticve form, we need to non-dimensionalize the primitive variables for removing the constraints of dimensions to analyse the numerical solution. It will certainly ease the method of analysis for a given problem.
Non-dimensional variables:
1. Temperature (T′T'):
T′=TT0T'=TT0
2. Density (ρ′ρ'):
ρ′=ρρ0ρ'=ρρ0
3. Length (x′x'):
x′=xLx'=xL
4. Speed of sound in the reservoir (a0a0):
a0=√γRT0a0=√γRT0
5. Velocity (V′V'):
V′=Va0V'=Va0
6. Time (t′t'):
t′=tLa0t'=tLa0
7. Area (A′A'):
A′=AA∗A'=AA∗
Conservative form:
1. Continuity equation:
∂(ρ′A′)∂t′+∂(ρ′A′V′)∂x′=0∂(ρ'A')∂t'+∂(ρ'A'V')∂x'=0
2. Momentum equation:
∂(ρ′A′V′)∂t′+∂[ρ′A′V′2+(1γ)ρ′A′]∂x′=1γ⋅p′⋅∂A′∂x′∂(ρ'A'V')∂t'+∂[ρ'A'V'2+(1γ)ρ'A']∂x'=1γ⋅p'⋅∂A'∂x'
3. Energy equation:
∂(ρ′(e′γ-1+(γ2)⋅V′2)A′)∂t′+∂(ρ′(e′γ-1+(γ2)⋅V′2)V′⋅A′+ρ′A′V′)∂x′=0∂(ρ'(e'γ−1+(γ2)⋅V'2)A')∂t'+∂(ρ'(e'γ−1+(γ2)⋅V'2)V'⋅A'+ρ'A'V')∂x'=0
The above is the generic form of non-dimensional conservation form of the continuity, momentum and energy equations for quasi-one-dimensional flow. Let us define the above equations in a simple generic form, which will be used in calculations instead of the above rigid and complex form. Elements of the solutions vector U, the flux vector F, and the source term J as follows,
U1=ρ′A′U1=ρ'A'
U2=ρ′A′V′U2=ρ'A'V'
U3=ρ′(e′γ-1+γ2V′2)A′U3=ρ'(e'γ−1+γ2V'2)A'
F1=ρ′A′V′F1=ρ'A'V'
F2=ρ′A′V′2+1γp′A′F2=ρ'A'V'2+1γp'A'
F3=ρ′(e′γ-1+γ2V′2)V′A′+p′A′V′F3=ρ'(e'γ−1+γ2V'2)V'A'+p'A'V'
J2=1γp′∂A′∂x′J2=1γp'∂A'∂x'
Before setting up numerical solution for conservative form, we must keep in mind that in the conservation form of the governing equations, the dependent variables (the variables for which we directly obtain numbers) are not primitive variables. In order to calculate the primitive variables, we need to decode the dependent variables and update our primitive variables in further steps.
Non-conservative form:
1. Continuity equation:
∂ρ′∂t′=-ρ′⋅∂V′∂x′-ρ′V′∂(lnA′)∂x′-V′∂ρ′∂x′∂ρ'∂t'=−ρ'⋅∂V'∂x'−ρ'V'∂(lnA')∂x'−V'∂ρ'∂x'
2. Momentum equation:
∂V′∂t′=-V′⋅∂V′∂x′-1γ⋅(∂T′∂x′+T′ρ′⋅∂ρ′∂x′)∂V'∂t'=−V'⋅∂V'∂x'−1γ⋅(∂T'∂x'+T'ρ'⋅∂ρ'∂x')
3. Energy equation:
∂T′∂t′=-V′⋅∂T′∂x′-(γ-1)⋅T′(∂V′∂x′+V′⋅∂(lnA′)∂x′)∂T'∂t'=−V'⋅∂T'∂x'−(γ−1)⋅T'(∂V'∂x'+V'⋅∂(lnA')∂x')
MacCormack Method:
MacCormack's method is an explicit finite-difference technique which is second order accurate in both space and time. MacCormack method generally uses three step approach(Predictor, Corrector and averaging) to calculate a final numerical solution for an iteration. In predictor step, forward difference is taken for spatial derivatives whereas in corrector step, rearward (backward) difference is taken for spatial derivatives. Furthermore, the solutions obtained from these two steps are averaged and final numerical solution is calculated for a single iteration.
Following is the step by step implementation of MacCormack method in non-conservative form of governing equations.
a) Predictor step:
1. Continuity equation:
(∂ρ∂t)ti=-ρti⋅Vti+1-Vti△x-ρtiVti⋅lnAi+1-lnAi△x-Vti⋅ρti+1-ρti△x→1(∂ρ∂t)ti=−ρti⋅Vti+1−Vti△x−ρtiVti⋅lnAi+1−lnAi△x−Vti⋅ρti+1−ρti△x→1
2. Momentum equation:
(∂V∂t)ti=-Vti⋅Vti+1-Vti△x-(1γ)⋅Tti+1-Tti△x+Ttiρti⋅ρti+1-ρti△x→2(∂V∂t)ti=−Vti⋅Vti+1−Vti△x−(1γ)⋅Tti+1−Tti△x+Ttiρti⋅ρti+1−ρti△x→2
3. Energy equation:
(∂T∂t)ti=-Vti⋅Tti+1-Tti△x-(γ-1)Tti⋅(Vti+1-Vti△x+Vti⋅lnAi+1-lnAi△x)→3(∂T∂t)ti=−Vti⋅Tti+1−Tti△x−(γ−1)Tti⋅(Vti+1−Vti△x+Vti⋅lnAi+1−lnAi△x)→3
By solving this equation, we obtain the predicted values of ρρ,VVandTT, denoted by the barred quantities,
ˉρt+△ti=ρti+(∂ρ∂t)ti△t→4¯ρt+△ti=ρti+(∂ρ∂t)ti△t→4
ˉVt+△ti=Vti+(∂V∂t)ti△t→5¯¯¯Vt+△ti=Vti+(∂V∂t)ti△t→5
ˉTt+△ti=Tti+(∂T∂t)ti△t→6¯¯¯Tt+△ti=Tti+(∂T∂t)ti△t→6
In equations 4 to 6, ρtiρti, VtiVtiand TtiTtiare known values at time t. And the values for the time derivatives in equation 4 to 6 are supplied directly by the equations 1 to 3.
Moving to the corrector step, taking backward difference in spatial derivatives and using predicted barred quantities
b) Corrector step:
1. Continuity equation:
(¯∂ρ∂t)t+△ti=-ˉρt+△ti⋅ˉVt+△ti-Vt+△ti-1△x-ˉρt+△ti⋅ˉVt+△ti⋅lnAi-lnAi-1△x-ˉVt+△ti⋅ˉρt+△ti-ˉρt+△ti-1△x→7(¯¯¯¯¯∂ρ∂t)t+△ti=−¯ρt+△ti⋅¯¯¯Vt+△ti−Vt+△ti−1△x−¯ρt+△ti⋅¯¯¯Vt+△ti⋅lnAi−lnAi−1△x−¯¯¯Vt+△ti⋅¯ρt+△ti−¯ρt+△ti−1△x→7
2. Momentum equation:
(¯∂V∂t)t+△ti=-ˉVt+△ti⋅ˉVt+△ti-ˉVt+△ti-1△x-(1γ)⋅ˉTt+△ti-ˉTt+△ti-1△x+ˉTt+△tiˉρt+△ti⋅ˉρt+△ti-ˉρt+△ti-1△x→8(¯¯¯¯¯¯∂V∂t)t+△ti=−¯¯¯Vt+△ti⋅¯¯¯Vt+△ti−¯¯¯Vt+△ti−1△x−(1γ)⋅¯¯¯Tt+△ti−¯¯¯Tt+△ti−1△x+¯¯¯Tt+△ti¯ρt+△ti⋅¯ρt+△ti−¯ρt+△ti−1△x→8
3. Energy equation:
(¯∂T∂t)t+△ti=-ˉVt+△ti⋅ˉTt+△ti-ˉTt+△ti-1△x-(γ-1)ˉTt+△ti⋅(ˉVt+△ti-ˉVt+△ti-1△x+ˉVt+△ti⋅lnAi-lnAi-1△x)→9(¯¯¯¯¯¯∂T∂t)t+△ti=−¯¯¯Vt+△ti⋅¯¯¯Tt+△ti−¯¯¯Tt+△ti−1△x−(γ−1)¯¯¯Tt+△ti⋅(¯¯¯Vt+△ti−¯¯¯Vt+△ti−1△x+¯¯¯Vt+△ti⋅lnAi−lnAi−1△x)→9
After solving the barred quantities from corrected step, we can calculate the average time derivatives,
c) Averaging time derivatives:
(∂ρ∂t)av=0.5⋅[(∂ρ∂t)ti+(¯∂ρ∂t)t+△ti]→10(∂ρ∂t)av=0.5⋅⎡⎣(∂ρ∂t)ti+(¯¯¯¯¯∂ρ∂t)t+△ti⎤⎦→10
(∂V∂t)av=0.5⋅[(∂V∂t)ti+(¯∂V∂t)t+△ti]→11(∂V∂t)av=0.5⋅[(∂V∂t)ti+(¯¯¯¯¯¯∂V∂t)t+△ti]→11
(∂T∂t)av=0.5⋅[(∂T∂t)ti+(¯∂T∂t)t+△ti]→12(∂T∂t)av=0.5⋅[(∂T∂t)ti+(¯¯¯¯¯¯∂T∂t)t+△ti]→12
Finaly, we have the corrected values of the flow-field variables at time (t+△tt+△t), we get,
ρt+△ti=ρti+(∂ρ∂t)av△t→13ρt+△ti=ρti+(∂ρ∂t)av△t→13
Vt+△ti=Vti+(∂V∂t)av△t→14Vt+△ti=Vti+(∂V∂t)av△t→14
Tt+△ti=Tti+(∂T∂t)av△t→15Tt+△ti=Tti+(∂T∂t)av△t→15
These variables are non-dimensional values and constitute our second eschelon of equations i.e the finite difference expressions of the governing equations in a form that pertains to MacCormack technique.
Note: Similarly, implement these steps for conservative form of the governing equations.
Initial and boundary conditions:
To start the time-marching calculations, we must stipulate initial conditions for ρρ, TT and VVas a function of x; that is, we must set up values of ρρ,TTand VVat time t = 0.
Non-conservative form:
Initial conditions:
ρ=1-0.3146xρ=1−0.3146x
T=1-0.2314xT=1−0.2314x
V=(0.1+1.09x)T12V=(0.1+1.09x)T12
Boundary conditions:
In our case, we are having a subsonic flow in the convergent section (M<1), choked flow at the throat section (M=1) and supersonic outflow in the divergent section (M>1). So by using method of characteristics for an unsteady, one-dimensional flow we can precisely allot inflow and outflow boundary conditions and can identify which quantities to state at the boundaries.
Fig (b) : Study of boundary conditions: subsonic inflow and supersonic outflow
What does these characteristics curve and streamline direction tells us?
- If any characteristic line is entering into the domain, implies that the value of one dependent flow field variable must be specified at that boundary.
- However, if any characteristic line is exiting the domain, implies that the value of the another dependent flow field variable must be allowed to float at that boundary i.e. it must be calculated in steps of time as function time wise solution of the flow field.
- Similarly, if streamline flows into the domain then specify third dependent flow field variable and if streamline flows out of the domain then third dependent flow field variable must be allowed to float.
In our case, at i = 1(inlet), we have primitive varibles like ρρ,AA,VV and TT. So which variables must be taken for initial and boundary conditions? Here, we know the ρρ and TT at node 1 are equal to the reservoir's ρρ and TT which is assumed to be ρ=1ρ=1 and `T = 1, so it becomes,
Inflow boundary condition:
ρ(1)=1ρ(1)=1 and T(1)=1T(1)=1 (Since it is Dirichlet boundary condition)
However, velocity and area tends to infinity, so for velocity, we need to calculate it from the information obtained from inner nodes,
V(1)=2⋅V(2)-V(3)V(1)=2⋅V(2)−V(3) (It is Neumann boundary condition)
Outflow boundary condition:
At the outflow boundary condition (Node =n), as we can see that, both the characteristics lines are out of the domain and streamline also flows out of the domain. This is because, at the outflow condition, the flow is supersonic (M>1), which is certainly greater than the speed of sound. So in this condition, information transport is very fast. Therefore, there are no flow field variables which require their values to be stipulated at the supersonic outflow boundary; all variables must be allowed to float at this boundary.
ρ(n)=2⋅ρ(n-1)-ρ(n-2)ρ(n)=2⋅ρ(n−1)−ρ(n−2)
V(n)=2⋅V(n-1)-V(n-2)V(n)=2⋅V(n−1)−V(n−2)
T(n)=2⋅T(n-1)-T(n-2)T(n)=2⋅T(n−1)−T(n−2)
The supersonic outflow boundary is purely Neumann boundary condition, which is extrapolating the boundary nodes using the information from inner nodes.
Conservative form:
In conservative form, we will initialize the nozzle profile rationally. As in conservative form, seen above from the governing equations, we have solution vectors, flux terms and a source term. So it is the prior matter of concern to state initial and boundary conditions in those terms rather than that of primitive variables.
Initial conditions:
These variations are slightly more realistic than those assumed in the nozzle shape and initial conditions in non-conservative form of governing equations discussed above. This is an 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 improved initial conditions than those discussed above in non-conservative form.
Assumptions for initial profile of the nozzle:
- For 0≤x′≤0.50≤x'≤0.5
ρ′=1ρ'=1
T′=1T'=1
- For 0.5≤x′≤1.50.5≤x'≤1.5
ρ′=1.0-0.366(x′-0.5)ρ'=1.0−0.366(x'−0.5)
T′=1.0-0.167(x′-0.5)T'=1.0−0.167(x'−0.5)
- For 1.5≤x′≤3.51.5≤x'≤3.5
ρ′=0.634-0.3879(x′-1.5)ρ'=0.634−0.3879(x'−1.5)
T′=0.833-0.3507(x′-1.5)T'=0.833−0.3507(x'−1.5)
Before stating the initial conditions in terms of solution vectors, we need to initialize velocity (V′V'). We have one of the dependent variable in our governing equation as U2U2, which is physically local mass flow; that is U2=ρ′A′V′U2=ρ'A'V'. Therefore, for the initial conditions only; let us assume a constant mass flow through the nozzle and calculate VV`as,
V′=U2ρ′A′=0.59ρ′A′V'=U2ρ'A'=0.59ρ'A'
Note: The value 0.59 is chosen for U2U2, because it is close to exact analytical value of the steady-state mass flow (which is for this case is 0.579)
Therefore, the initial condition for V′V'as a function of x′x' is obtained by substituing the ρ′ρ'and A′A' variation. Finally, the initial conditions for U1U1, U2U2, U3U3 are obtained by substituting the above variations for ρ′ρ', T′T' and V′V'into the definitions discusssed above.
U1=ρ′A′U1=ρ'A'
U2=ρ′A′V′U2=ρ'A'V'
U3=ρ′(e′γ-1+γ2V′2)A′U3=ρ'(e'γ−1+γ2V'2)A'
where, e′=T′e'=T'. Certainly, for the initial conditions described above, V′V'is calculated such that U2=ρ′A′V′=0.59U2=ρ'A'V'=0.59
Boundary conditions:
In conservative form, by referring fig (b) and absorbing the same method of characteristics approach we can allot boundary conditions for inflow and outflow boundary conditions. However, in this case, the boundary conditions will be in terms of solution vectors. Similarly as discussed above, at the subsonic inflow (M<1), we need to specify two quantities at the initial node, where we have ρ′ρ',A′A' and U2=ρ′A′V′=0.59U2=ρ'A'V'=0.59(Mass flow rate) specified and one quantity will be allowed to float by referring left hand characteristic line. However, at the supersonic outflow (M>1), the information is travelling exceeds the speed of sound, hence it is irrational to specify outflow boundary condition. Therefore, all the dependent variables are allowed to float and they will be calculated in steps of time as a function time-wise solution of the flow field.
Inflow boundary condition:
The inflow boundary condition can be given as,
U1=ρ′(1)A′(1)U1=ρ'(1)A'(1) (Dirichlet boundary condition)
U2=2⋅U2(2)-U2(3)U2=2⋅U2(2)−U2(3)
V(1)=U2(1)U1(1)V(1)=U2(1)U1(1)
Note: This velocity at the first node will be substituted below in calculating U3U3
U3=U1(1)[T(1)γ-1+γ2(V′(1))2]U3=U1(1)[T(1)γ−1+γ2(V'(1))2] (Neumann boundary condition)
Outflow boundary condition:
All the dependent variables will be allowed to float. Hence, the outflow boundary can be given as,
U1(n)=2⋅U1(n-1)-U1(n-2)U1(n)=2⋅U1(n−1)−U1(n−2) (Neumann boundary condition)
U2(n)=2⋅U2(n-1)-U2(n-2)U2(n)=2⋅U2(n−1)−U2(n−2) (Neumann boundary condition)
U3(n)=2⋅U3(n-1)-U3(n-2)U3(n)=2⋅U3(n−1)−U3(n−2) (Neumann boundary condition)
Time step calculation:
The governing system of equations is hyperbolic with respect to time. So there is stability constraint on this system which is known as Courant-Friedrich-Lewy (CFL) condition. The Courant number is a simple stability analysis of a linear hyperbolic equations that gives result C≤1C≤1 for an explicit numerical solution to be stable.
CFL condition can be given as,
CFL condition = (any quantity)*(time-step with which information would travel)/(element step size)
So, in our case, the quantity which is going to travel is velocity, or in other sense velocity is the quantity which will be travelling with certain time from one node to other. So here velocity of the fluid at a point in the flow is coupled with local speed of sound. The quantity can be expressed to be as a+V′a+V'
CFL = a + V'⋅△t△x⋅△t△x
We will assume the Courant number as C = 0.5, accordingly will calculate the △t. In addition to it, we will calculate (△t)ti
at all the grid points, i = 1 to i = n, and then we'll choose the minimum value of △t in our calculations.
After considering all these governing equations and their respective pre-requisites to solve the problem, then we will move towards to the MATLAB code for execution.
Main MATLAB code:
clear all
close all
clc
% Defining inputs
L = 3; % Length of the domain
n = 31; % Total number of nodes in x-dimension
x = linspace(0,L,n); % X-dimension grids
dx = x(2) - x(1); % Elemental step size
gamma = 1.4; % Specific heat capacity of air
a = 1+2.2*(x-1.5).^2; % Area of the nozzle
cfl = 0.5; % CFL Number
nt = 1400; % Total time steps
time_step = linspace(1,nt,nt); % Dummy variable of total number of time steps
% Calling conservative function
[rho_1,rho_th_1,t_1,t_th_1,v_1,p_1,p_th_1,mach_1,mach_th_1,m_fl_1,k_1,m_fl_an_1] = conserv(x,n,dx,gamma,a,cfl,nt);
% Calling non-conservative function
[rho_2,rho_th_2,t_2,t_th_2,v_2,p_2,p_th_2,mach_2,mach_th_2,m_fl_2,k_2,m_fl_an_2] = non_conserv(x,n,dx,gamma,a,cfl,nt);
% Plotting results
% Comparison of timewise variation of flow field variables at the throat
% section
figure(3)
subplot(4,1,1)
plot(time_step,rho_th_1)
hold on
plot(time_step,rho_th_2)
xlabel('Number of time steps')
ylabel('Density (rho)')
legend('Non-conservation form','Conservation form')
title('Time-wise variation of flow variables at throat s/c (Density)')
subplot(4,1,2)
plot(time_step,t_th_1)
hold on
plot(time_step,t_th_2)
xlabel('Number of time steps')
ylabel('Temperature (t)')
legend('Non-conservation form','Conservation form')
title('Time-wise variation of flow variables at throat s/c (Temperature)')
subplot(4,1,3)
plot(time_step,p_th_1)
hold on
plot(time_step,p_th_2)
xlabel('Number of time steps')
ylabel('Pressure (p)')
legend('Non-conservation form','Conservation form')
title('Time-wise variation of flow variables at throat s/c (Pressure)')
subplot(4,1,4)
plot(time_step,mach_th_1)
hold on
plot(time_step,mach_th_2)
xlabel('Number of time steps')
ylabel('Mach number (M)')
legend('Non-conservation form','Conservation form')
title('Time-wise variation of flow variables at throat s/c (Mach Number)')
% Plots of variation of flow paramters inside the nozzle
figure(4)
subplot(4,1,1)
plot(x,rho_1)
hold on
plot(x,rho_2)
xlabel('Non-dimensional distance through the nozzle(x)')
ylabel('Density (rho)')
legend('Non-conservation form','Conservation form')
title('Non dimensional density variation (rho) along x at differrent time steps')
subplot(4,1,2)
plot(x,t_1)
hold on
plot(x,t_2)
xlabel('Non-dimensional distance through the nozzle(x)')
ylabel('Temperature (t)')
legend('Non-conservation form','Conservation form')
title('Non dimensional temperature variation (t) along x at differrent time steps')
subplot(4,1,3)
plot(x,p_1)
hold on
plot(x,p_2)
xlabel('Non-dimensional distance through the nozzle(x)')
ylabel('Pressure(p)')
legend('Non-conservation form','Conservation form')
title('Non dimensional pressure variation (p) along x at differrent time steps')
subplot(4,1,4)
plot(x,mach_1)
hold on
plot(x,mach_2)
xlabel('Non-dimensional distance through the nozzle(x)')
ylabel('Mach number (M)')
legend('Non-conservation form','Conservation form')
title('Non dimensional Mach number variation (M) along x at differrent time steps')
% Plot steady state mass flow variations obtained with non-conservation and
% conservation form of the governing equations
figure(5)
plot(x,m_fl_1,'r')
hold on
plot(x,m_fl_2,'b')
plot(x,m_fl_an_1,'g')
xlabel('Non-dimensional distance through the nozzle(x)')
ylabel('Non dimensional mass flow rate (rho*v*a)')
legend('Conservative form',' Non-conservative form','Exact analytical solution')
title('Comparison of mass flow rate variation')
% To calculate the performance of the functions in terms of time required
% for execution.
% Conservative form
myfunc_con = @()conserv(x,n,dx,gamma,a,cfl,nt)
t_con= timeit(myfunc_con);
% Non conservative form
myfunc_noncon = @() non_conserv(x,n,dx,gamma,a,cfl,nt)
t_noncon= timeit(myfunc_noncon);
% Displaying time required to execute the function
fprintf('Time taken to execute Conservative form = %f',t_con)
fprintf('n Time taken to execute Non-conservative form = %f',t_noncon)
Function code for Conservative form:
function [rho,rho_th,t,t_th,v,p,p_th,mach,mach_th,m_fl,k,m_fl_an] = conserv(x,n,dx,gamma,a,cfl,nt)
% Initialize variables
for s = 1:n
if x(s)>=0 && x(s)<=0.5
rho(s)= 1;
t(s)= 1;
else
if x(s)>=0.5 && x(s)<=1.5
rho(s) = 1-0.366*(x(s)-0.5);
t(s) = 1 - 0.167*(x(s)-0.5);
else
if x(s)>=1.5 && x(s)<=3.5
rho(s) = 0.634-0.3879*(x(s)-1.5);
t(s) = 0.833 - 0.3507*(x(s)-1.5);
end
end
end
end
% Assume constant mass flow through the nozzle
m_fl = 0.590; % Constant mass flow rate
v = m_fl./(rho.*a); % Velocity at constant mass flow rate
% To calculate minimum time step using CFL condition
for g = 1:n
del_t(g) = cfl*(dx/(t(g)^0.5 + v(g)));
end
dt = min(del_t);
% Time step used in John.D.Anderson is (dt = 0.0267). However, this dt is
% not achieved using CFL loop. Therefore extracting minimum value of dt
% for the calculation (dt = 0.0201)
% Define solution vector (u)
u1 = rho.*a;
u2 = rho.*a.*v;
u3 = u1.*(t./(gamma-1) + ((gamma/2).*(v.^2)));
% Time loop
for k = 1:nt
% Initialize the primitive variables
u1_old = u1;
u2_old = u2;
u3_old = u3;
% Define flux vector(f)
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));
% Steady state analytical mass flow distribution
m_fl_an = 0.579*ones(n);
%% Predictor step
for i = 2:n-1
% Predefined variables
df1dx = (f1(i+1)-f1(i))/dx;
df2dx = (f2(i+1)-f2(i))/dx;
df3dx = (f3(i+1) - f3(i))/dx;
dadx = (a(i+1) - a(i))/dx;
j2_p = (1/gamma).*rho(i).*t(i).*dadx;
% Continuity equation
du1dt_p (i) = -df1dx;
% Momentum equation
du2dt_p(i) = -df2dx + j2_p;
% Energy equation
du3dt_p(i) = -df3dx;
% Solution update
u1(i) = u1(i) + du1dt_p(i)*dt;
u2(i) = u2(i) + du2dt_p(i)*dt;
u3(i) = u3(i) + du3dt_p(i)*dt;
end
% Primitive variables update
rho = u1./a;
v = u2./u1;
t = (gamma-1)*((u3./u1)-((gamma/2)*((u2./u1).^2)));
%Decode flux variables for corrector step
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 i = 2:n-1
% Predefined variables
df1dx = (f1(i)-f1(i-1))/dx;
df2dx = (f2(i)-f2(i-1))/dx;
df3dx = (f3(i) - f3(i-1))/dx;
dadx = (a(i) - a(i-1))/dx;
j2_c = (1/gamma).*rho(i).*t(i).*dadx;
% Continuity equation
du1dt_c (i) = -df1dx;
% Momentum equation
du2dt_c(i) = -df2dx + j2_c;
% Energy equation
du3dt_c(i) = -df3dx;
end
% Compute the average time derivative
du1dt = 0.5*(du1dt_p + du1dt_c);
du2dt = 0.5*(du2dt_p + du2dt_c);
du3dt = 0.5*(du3dt_p + du3dt_c);
% Final solution update
for h = 2:n-1
u1(h) = u1_old(h) + du1dt(h)*dt;
u2(h) = u2_old(h) + du2dt(h)*dt;
u3(h) = u3_old(h) + du3dt(h)*dt;
end
% Calculating primitive variables with the help of solution vector
for o = 2:n-1
rho(o) = u1(o)/a(o);
v(o) = u2(o)/u1(o);
t(o) = (gamma-1)*(u3(o)/u1(o) - (gamma/2)*v(o)^2);
end
% Applying boundary conditions
% Inlet
u1(1) = rho(1)*a(1); % Dirichlet BC
u2(1) = 2*u2(2) - u2(3); % Neumann BC
v(1) = u2(1)/u1(1);
u3(1) = u1(1)*((t(1))/(gamma-1) + (gamma/2)*v(1)^2); % Neumann BC
% Outlet
u1(n) = 2*u1(n-1) - u1(n-2); %Neumann BC
u2(n) = 2*u2(n-1) - u2(n-2); %Neumann BC
u3(n) = 2*u3(n-1) - u3(n-2); %Neumann BC
% Calculating other variables
m_fl = rho.*a.*v; % Mass flow rate
mach = v./((t).^0.5); % Mach number
p = rho.*t; % Pressure
% Calculating throat variables
rho_th (k) = rho(15);
t_th (k) = t(15);
p_th (k) = p(15);
mach_th (k) = mach(15);
% Plotting time-wise variation of mass flow rate
if k == 50
figure(1)
plot(x,m_fl,'r','linewidth',2)
hold on
else
if k == 200
plot(x,m_fl,'b','linewidth',2)
else
if k == 400
plot(x,m_fl,'m','linewidth',2)
else
if k == 800
plot(x,m_fl,'g','linewidth',2)
else
if k == 1400
plot(x,m_fl,'y','linewidth',2)
plot (x,m_fl_an,':','linewidth',2)
hold off
end
end
end
end
end
grid on
xlabel('Non-dimensional distance through the nozzle(x)')
ylabel('Non-dimensional mass flow rate (mfl)')
end
title('Non dimensional mass flow rate (Conservative form) along x at differrent time steps')
legend('50 time-step','200 time-step','400 time-step','800 time-step','1400 time-step','Exact analytical solution','location','northwest')
end
Function code for Non-conservative form:
function [rho,rho_th,t,t_th,v,p,p_th,mach,mach_th,m_fl,k,m_fl_an] = non_conserv(x,n,dx,gamma,a,cfl,nt)
% Initial condition
rho = 1- (0.3146*x); % Density
t = 1 - (0.2314*x); % Temperature
v = (0.1 + 1.09*x).*(t).^0.5; % Velocity
% To calculate minimum time step using CFL condition
for g = 1:n
del_t(g) = cfl*(dx/(t(g)^0.5 + v(g)));
end
dt = min(del_t);
% Time step used in John.D.Anderson is (dt = 0.0267). However, this dt is
% not achieved using CFL loop. Therefore extracting minimum value of dt
% for the calculation (dt = 0.0201)
% Time loop
for k = 1:nt
% Initialize primitive variables
rho_old = rho;
v_old = v;
t_old = t;
% Steady state analytical mass flow distribution
m_fl_an = 0.579*ones(n);
%% Predictor step
for i = 2:n-1
% Predefined variables
drhodx = (rho(i+1)-rho(i))/dx ;
dvdx = (v(i+1) - v(i))/dx;
dlogadx = (log(a(i+1)) - log(a(i)))/dx;
dtdx = (t(i+1) - t(i))/dx;
% Continuity equation
drhodt_p (i) = -(rho(i)*dvdx) - (rho(i)*v(i)*dlogadx) - v(i)*(drhodx);
% Momentum equation
dvdt_p(i) = -v(i)*dvdx - (1/gamma)*(dtdx+(t(i)/rho(i))*drhodx);
% Energy equation
dtdt_p(i) = -(v(i)*dtdx) - (gamma-1)*t(i)*(dvdx + v(i)*dlogadx);
% Solution update
rho(i) = rho(i) + drhodt_p(i)*dt;
v(i) = v(i) + dvdt_p(i)*dt;
t(i) = t(i) + dtdt_p(i)*dt;
end
%% Corrector step
for m = 2:n-1
% Predefined variables
drhodx = (rho(m)-rho(m-1))/dx ;
dvdx = (v(m) - v(m-1))/dx;
dlogadx = (log(a(m)) - log(a(m-1)))/dx;
dtdx = (t(m) - t(m-1))/dx;
% Continuity equation
drhodt_c (m) = -(rho(m)*dvdx) - (rho(m)*v(m)*dlogadx) - v(m)*(drhodx);
% Momentum equation
dvdt_c(m) = -v(m)*dvdx - (1/gamma)*(dtdx+(t(m)/rho(m))*drhodx);
% Energy equation
dtdt_c(m) = -(v(m)*dtdx) - (gamma-1)*t(m)*(dvdx + v(m)*dlogadx);
end
%% Compute the average time derivative
drhodt = 0.5*(drhodt_p + drhodt_c);
dvdt = 0.5*(dvdt_p + dvdt_c);
dtdt = 0.5*(dtdt_p + dtdt_c);
% Final solution update
for l = 2:n-1
rho(l) = rho_old(l) + drhodt(l)*dt;
v(l) = v_old(l) + dvdt(l)*dt;
t(l) = t_old(l) + dtdt(l)*dt;
% To calculate mass flow rate
m_fl(l) = rho(l)*a(l)*v(l);
end
%% Applying boundary conditions
% Inlet
rho(1)= 1; % Dirichlet BC
t(1) = 1; % Dirichlet BC
v(1) = 2*v(2) - v(3); % Neumann BC
% Outlet
v(n) = 2*v(n-1) - v(n-2); % Neumann BC
rho(n) = 2*rho(n-1) - rho(n-2); % Neumann BC
t(n) = 2*t(n-1) - t(n-2); % Neumann BC
% Calculating other variables
m_fl = rho.*a.*v; % Mass flow rate
mach = v./((t).^0.5); % Mach number
p = rho.*t; % Pressure
% Calculating throat variables
rho_th (k) = rho(15);
t_th (k) = t(15);
p_th (k) = p(15);
mach_th (k) = mach(15);
% Plotting time-wise variation of mass flow rate
if k == 50
figure(2)
plot(x,m_fl,'r','linewidth',2)
hold on
else
if k == 200
plot(x,m_fl,'b','linewidth',2)
else
if k == 400
plot(x,m_fl,'m','linewidth',2)
else
if k == 800
plot(x,m_fl,'g','linewidth',2)
else
if k == 1400
plot(x,m_fl,'y','linewidth',2)
plot(x,m_fl_an,':','linewidth',2)
hold off
end
end
end
end
end
grid on
xlabel('Non-dimensional distance through the nozzle(x)')
ylabel('Non-dimensional mass flow rate (mfl)')
end
title('Non dimensional mass flow rate (Non-conservative form) along x at differrent time steps')
legend('50 time-step','200 time-step','400 time-step','800 time-step','1400 time-step','Exact analytical solution','location','northwest')
end
Results and discussions:
Plots:
1. Steady state distribution of primitive variables inside the nozzle:
At nodes = n = 31:
At nodes = n = 61:
Discussion:
The steady state distribution inside the nozzle is resulted as above. After 1400 time steps, the numerically calculated primitive variables are converged to the exact analytical solution. The variations along the non-dimensional distance of the nozzle is explicitly resulted. Initially some of the primitive variables were set to initial condition and other were allowed to float as a part of inflow subsonic boundary condition. At n = 1, non-dimensional pressure or pressure ratio (pp0) of the fluid element was set to 1 (reservoir end), then as a part of numerical explicit scheme (MacCormack method), the pressure variations were calculated as function of time along the nozzle. There is decrease in the pressure as the flow changes from subsonic (convergent s/c, M<1) to sonic (throat s/c, M=1) state and similarly from sonic (throat s/c, M=1) to supersonic (divergent s/c, M>1) state. Similar is the case for density and temperature. The pressure difference present inside the system and in the surronding is responsible for the fluid flow. This pressure difference (reservoir and exit/surronding) gradually decreases as a function of time and length of the nozzle as the fluid travels towards the exit of the divergent section. Conversely, velocity of the fluid will increase travelling through convergent (subsonic, a-V), throat(sonic, V = a) and divergent section(supersonic, V+a). Similarly, as the velocity increases, Mach number will increase accordingly as function of length of the nozzle. Both conservative and non-conservative form of the governing equations show steady state results and show similar results at the macroscopic level. Also there are similar results obtained using grid points (nodes) as 31 and 61.
2. Time-wise variations of the primitive variables:
At nodes = n = 31:
At nodes = n = 61:
Discussion:
The variation of all the primitive variables at the throat section are resulted as function of time. The results are changing continuously as the time advances from an initial state (t=0). The MacCormack method is converging numerical results with each iteration in fast and robust manner. As we know, the throat section is responsible for increasing the velocity of the fluid element to a sonic state. It is prime reference for our calculations to compare results at the throat section with exact analytical solution at throat. So it neccessary to plot variations of the primitive variables as a part of time marching process. When using nodes as 31, there is very minute difference between variations of the primitive variables at the end of time (1400 time-steps) using conservative and non-conservative form. However, when the grid size is doubled, both the conservative and non-conservative form calculating the primitive variables clearly converges the numerical solution at approx 1100th time step. However, there is very slight and neglible difference in between both the conditions as both the curves obtain steady state at the end of 1400 time steps (i.e. numerically converged to exact analytical solution).
3. Variation of mass flow distribution inside the nozzle at different time-steps during the time-marching process:
3.1.1 Conservative form using nodes = n = 31:
3.2.1 Conservative form using Nodes = n = 61:
3.1.2 Non-conservative form using nodes = n = 31:
3.2.2 Non-conservative form Nodes = n = 61:
Discussion:
In these above plots, the mass flow distribution has resulted as function of advancement in time-step. Irrespective of the form of governing equation used, the mass flow distribution at first 50 time-step shows drastic variation. As the time step progresses, the mass flow distribution variation decreases and somewhat linearizes eventually. This is because the numerical result takes time to converge to the exact analytical solution. In addition to it, the conservative form has less variations among the curves as compared to non-conservative form. This is because we provided improved initial condition for conservative form. Due to this, the initial condition was defined close to steady state solution which refined the process flow. Moreover, the results obtained from the solution with nodes (n) = 31 and nodes (n) = 61 are similar and show slight negligible differences in their numerical solution at the macroscopic level. Upon some magnification, certainly the final numerical solution with 61 nodes will show more closer results than that with 31 nodes. However, the difference between them is very less. So at 1400th time step, both the form of numerical solution (Conservative and non-conservative) converged to the exact analytical solution numerically.
4. Comparison of normalized mass flow rate distributions of both forms:
At nodes = n = 31:
At nodes = n = 61:
Discussion:
The comparison of mass flow distribution is plotted against the non-dimensional nozzle distance. By observing the above plot, we can state that the mass flow distribution predicted by the conservative form of equations shows good results than that of the non-conservative form. This is because of some inferences from the resulted graph. The conservation form gives a distribution that is much closer to being a constant. However, the non-conservative form shows significant variation and oscillations at the inflow and outflow boundaries. The reason behind conservation form to yield better results is that mass flow is one of the dependent variables in the equations, which mean that the mass flow is a primary result of the equations. In contrast to it, the dependent variables in the nonconservation form of the equations are the primitive variables, and the mass flow is obtained only as a secondary result.
Observation table:
Table 1: Steady-state results, using conservative form (n = 31):
Tabular indications:
Upper row in green colour = Variables at the inflow boundaries (Reservoir end)
Row in blue color = Variables at throat section
Lower row in green colour = Variables at the outflow boundaries (Exit of the divergent section)
Table 2: Steady-state results, using non-conservative form (n = 31):
Tabular indications:
Upper row in green colour = Variables at the inflow boundaries (Reservoir end)
Row in blue color = Variables at throat section
Lower row in green colour = Variables at the outflow boundaries (Exit of the divergent section)
Table 3: Comparison of numerical flow field variables with exact analytical solution:
Table 4: Estimation of simulation time:
Discussion:
The tabular representation is also significant because it gives the numeric value of any quantity at a particular node position. In table 1, we have all the values of primitive and dependent variables calculated at each node position using conservation form of equations. Similarly, in table 2 we have all the primitive variables at each node position using non-conservative form. By observing the throat values and comparing it with exact analytical values, we can evaluate both the form of equations. By referring table 3, it is clear that non-conservation form offers closest values of the primitive variables (pressure, temperature, density and mach number) to the exact analytical solution. Comparatively, conservation form of equations lacks in calculating primitive variables closest to the exact one. This might be because, in conservation form, the dependent variables are calculated first and then decoded to calculate the primitive variables which affects the accuracy of the numerical solution. However, in the non-conservation form of equations, dependent variables are the primitive variables which results into an accurate numerical solution for the primitive variables.
While observing results for different number of grids, one interesting aspect in numerical analysis was resulted that though the number of grid points are doubled there is no or slight change in the results. However, we know that, by increasing the number of grid points, we obtain accurate results than that for coarser grid. In our case, there is very slight variation (negligible) in our results though we doubled the grid points. This is because of the grid independence. In simple terms, grid independence means calculational results change so little along with a denser or looser grid that the truncation error can be ignored in numerical simulation.
The time required to produce/simulate the results is also an important aspect in CFD analysis. However, in our case, certainly for 61 grid points the time required is more than that for 31 grids. The differences are very negligible and can be ignored for this case. Again upon further magnification into the results we can state that the non-conservative form of equations consume less amount of time to produce/simulate the results as compared to the conservative form. This is because of the algorithm used in solving primitive variables in non-conservative form of equations is straightforward and simple as compared to conservative form of equations.
Conclusion:
In summary, for the given flow problem, we cannot establish an assertion that which form of equation (conservation or non-conservation) has superiority over the other form. In essence, from our previous discussions, we can only make the following observations:
a) The conservation form yields a better mass flow distribution. The conservation form does better job in conserving mass.
b) The non-conservation form of equations calculates more accurate primitive variables than that of the conservation form of equations. This is because of the straightforward numerical algorithm of non-conservation form. Which certainly relates to the residuals incurred in the numerical solution is comparatively less. This is often referred to as an index of "quality" of the numerical algorithm.
c) The conservation form of equations seems to produce slightly more accurate results for the flux variables such as U2=ρ′A′V′= Mass flow rate than that of the non-conservation form of equations.
d) The initial condition provided to conservation form of equations is slightly improved because of sensitive nature of the numerical algorithm, so it is certain to yield better results in some aspects of the numerical solution than that of the non-conservative form.
e) If we note in terms of calculational effort and complexity of computing to achieve a numerical solution, the conservation form requires comparatively more work. It is because of the fact that, we need to code dependent variables in terms of primitive variables, then decode the primitive variables from the flux variables; such decoding is not necessary for the non-conservation form.
f)The results in either case are certainly satisfactory. After reviewing all those thick and thin of the forms of governing equations, we cannot establish an assertion pertaining to the superiority of one form over the other.
References:
1. Computational fluid dynamics: The basics with applications by John D. Anderson, Jr.
2. https://www.grc.nasa.gov/www/k-12/airplane/nozzled.html
Leave a comment
Thanks for choosing to leave a comment. Please keep in mind that all the comments are moderated as per our comment policy, and your email will not be published for privacy reasons. Please leave a personal & meaningful conversation.
Other comments...
Be the first to add a comment
Read more Projects by Tanmay Pathak (15)
Week 10 - Simulating Combustion of Natural Gas.
Aim: To perform and study combustion simulation of a non-premixed model. Introduction and Theory: Combustion or burning is a high-temperature exothermic redox chemical reaction between a fuel(the reductant) and an oxidant, usually, atmospheric oxygen. Complete combustion is stoichiometric concerning the…
05 Aug 2021 12:02 PM IST
Numerical Analysis of Steady and Transient Two-Dimensional Heat Conduction Problem using MATLAB.
Aim: To solve steady and unsteady two dimensional heat conduction equation using implicit and explicit methods. Objectives: 1. Write a MATLAB code to solve steady state two dimensional heat conduction equation using implicit method. 2. Write a MATLAB code to solve unsteady (transient) state two dimensional…
14 Jul 2021 10:52 AM IST
Analysis and Validation of Quasi One Dimensional De Laval Nozzle using MacCormack Method.
Aim: Simulation of quasi 1D subsonic-supersonic nozzle flow using MacCormack method. Objectives: The objective of this project is to simulate an isentropic flow through a quasi 1D subsonic-supersonic nozzle using MacCormack method for conservative and non-conservative form of governing equations. Theory:…
14 Jul 2021 09:55 AM IST
Simulation of Flow Through a Pipe in OPENFOAM and Validation with the Analytical Results.
Aim: To simulate flow through a pipe in OPENFOAM and validate the results. Objectives: 1. To validate of hydrodynamic length with the numerical result. 2. Validate the fully developed flow velocity profile with its analytical profile. 3. Validate maximum velocity and pressure drop for fully developed…
14 Jul 2021 09:33 AM IST
Related Courses
Skill-Lync offers industry relevant advanced engineering courses for engineering students by partnering with industry experts.
Our Company
4th Floor, BLOCK-B, Velachery - Tambaram Main Rd, Ram Nagar South, Madipakkam, Chennai, Tamil Nadu 600042.
Top Individual Courses
Top PG Programs
Skill-Lync Plus
Trending Blogs
© 2025 Skill-Lync Inc. All Rights Reserved.