All Courses
All Courses
Courses by Software
Courses by Semester
Courses by Domain
Tool-focused Courses
Machine learning
POPULAR COURSES
Success Stories
OBJECTIVE: In this challenge, you are going to simulate the isentropic flow through a quasi 1D subsonic-supersonic nozzle. You will derive both the conservation and non-conservation forms of the governing equations and solve them using MacCormack's technique. You need to determine the steady-state temperature distribution…
Sagar Biswas
updated on 22 Sep 2022
OBJECTIVE: In this challenge, you are going to simulate the isentropic flow through a quasi 1D subsonic-supersonic nozzle. You will derive both the conservation and non-conservation forms of the governing equations and solve them using MacCormack's technique. You need to determine the steady-state temperature distribution for the flow-field variables and investigate the difference between the two forms of governing equations by comparing their solutions.
You need to show the following plots inside your report
1. Steady-state distribution of primitive variables inside the nozzle
2. Time-wise variation of the primitive variables
3. Variation of Mass flow rate distribution inside the nozzle at different time steps during the time-marching process
4. Comparison of Normalized mass flow rate distributions of both forms
Expected Results:
THEORY:
Conservative & Non-Conservative forms of Governing Equations:
Before getting into the theory for Macormack Method, we need to understand the concept and working principle of conservative and non-conservative forms of governing equations.
Both from the theoretical & mathematical perspectives, they're the same in terms of behavior when they're being solved numerically. Both forms suitably represent the fundamental principles of Mass conservation, Newton's 2nd law, and Energy conservation. By using Calculus, both forms can be converted into one another but when viewed from a numerical perspective, the choice or selection of one form over another would rather result in making a world of a difference.
Each form consists of 3 different equations representing the conservation of mass, energy & momentum. To better understand this, we're now considering the continuity equation derived for an infinitesimal control volume as shown below:
Non-conservative form: DρDt+ρ∇V=0DρDt+ρ∇V=0
Conservative form: ∂ρ∂t+∇(ρV)=0∂ρ∂t+∇(ρV)=0
Now, we're going to study the differences between these two forms of equations by comparing them with one another.
The first difference we can notice between these two forms is the 'Time derivative' term. The local derivative in the conservation form is replaced with a substantial derivative in the non-conservation. In the conservation form, the derivation is made under the assumption that the control volume is fixed and hence the variation of density with respect to time is provided by the local derivative while in the case of the non-conservation form, the derivation is made under the assumption that the control volume is moving. Here, the variation of density with respect to time is given by the sum of the local derivative(variation due to unsteady nature) and the convective derivative(variation due to moving control volume).
The 2nd difference that can be noticed is the position of the ∇∇ operator. It is called the differential/divergence operator whose dot product with a vector results in providing the divergence(volume intensity of outward flux) of that vector field. In the non-conservation form, it can be seen that the divergence operator is applied only to the velocity field and not density as density lies outside the operator while in conservation form, the divergence operator is applied on both density and velocity field and there is no variable outside the divergence operator.
In the continuity equation of the conservation form, the quantity on which the divergence operator is acting in mass flux is the rate of mass flow per unit area or the momentum per unit volume which makes sense only when the control volume is not moving.
There are two approaches that can be used to derive the governing equations in fluid dynamics. Those are Eulerian & Lagrangian approaches. In the Eulerian approach, we consider a fixed point in space and we record the properties of fluid elements passing through it while in the case of the Lagrangian approach, we follow the fluid particles as they travel through the flow.
Each form has its own features that make it suitable for some problems more than the others. In terms of programming convenience, the non-conservation form is better as it is easier to program because the equations are solved in terms of the primitive variables(ρρ, V, T) & their values can be obtained directly. In conservation form, the equations are expressed in flux variables and those are needed to be modified in terms of solution vectors. The equations are solved for the solution vectors and then the primitive variables are determined from them. This adds a fair amount of complexity while solving the governing equations in the conservation form.
The non-conservation form is considered to be more accurate between the two equations as the primitive variables are directly solved in the governing equations. In the conservation form, the primitive variables are only computed indirectly from the flux variables since the flux variables are directly solved for in the conservation equations which causes a loss in the order of accuracy of the primitive variables in the conservation form. It is important to note that an increased number of calculations will lead to more errors due to the computer rounding-off decimal values to machine precision. In that respect, if we consider the flux variables, however, the conservation equations are more accurate for the same reasons mentioned above.
The primitive variables converge to a numerical steady-state faster when using the non-conservation form because of the fact that the dependent variables are primitive in non-conservation forms but not in conservation form. Consequently, the simulation time needed for the non-conservation form is shorted than the conservation form.
The most important difference between the two forms comes in regard to their stability. We know that the equations we solve in CFD may not always be continuous. There might exist certain discontinuities in some regions of the domain due to the physical or mathematical nature of the equations we're solving. The primitive variables like density, temperature & velocity might experience a sudden jump across these discontinuities. When the governing equations are expressed using the primitive variables (like in non-conservation form), they may develop large gradients in the solution that might cause the solution to blow up. On the other hand, the variables possess an interesting property that they remain constant across shocks or discontinuities.
Here, we're considering a steady-state, isentropic flow through a convergent-divergent nozzle. The flow at the inlet to the nozzle comes from a reservoir where the pressure and temperature are denoted by P0P0and T0T0, respectively. The cross-sectional area of the reservoir is large (theoretically, A=∞A=∞), and hence the velocity is very small(V=0). Thus, P0P0and T0T0 are the stagnation values, or total pressure and total temperature, respectively. The stagnation properties at a point are defined as those which are to be obtained if the local flow were imagined to cease to zero velocity isentropically. The flow expands isentropically to supersonic speeds at the nozzle exit, where the exit pressure, temperature, velocity, and Mach number are denoted by PePeTeTe, VeVe & MeMe, respectively. The flow is locally subsonic in the convergent section of the nozzle, sonic at the throat(minimum area), and supersonic at the divergent section. The sonic flow(M=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*. Similarly, the sonic flow values of pressure and temperature are denoted by P* and T* respectively. The area of the sonic throat is denoted by A*. We assume that at a given section, where the cross-sectional area is A, the flow properties are uniform across that section. Hence, although the area of the nozzle changes as a function of distance along with the nozzle, x, and therefore in reality the flow field is two-dimensional (the flow varies in the two dimensional XY space), we make the assumption that the flow properties vary only with x; this is tantamount to assuming uniform flow properties across any given cross-section. Such flow is defined as the quasi-one-dimensional flow.
The governing continuity, momentum, and energy equations for this quasi-one-dimensional, steady, isentropic flow can be expressed, respectively, as:
Here, subscripts 1 and 2 denote different locations along with the nozzle.
The Mach number variation through the nozzle is governed exclusively by the area ratio A/A* through the relation:
Where, γγ= ratio of specific heats = CpCvCpCv. For standard air conditions, γγ= 1.4.
For a nozzle where A is specified as a function of x, hence A/A* is known as a function of x, then the above equation allows the implicit calculation of M as a function of x. The variation of pressure, density, and temperature as a function of Mach number (and hence as a function of A/A*, thus x) is given as,
Here, in this case, the force exerted on the gas to accelerate it through the nozzle is supplied by the pressure ratio across the nozzle, P0P0/PePe.
Also, when they're solved analytically, we arrive at the following results shown in the graphs below:
Figure: Qualitative aspects of quasi-one-dimensional nozzle flow: isentropic subsonic-supersonic solution.
For a nozzle with a specified area ratio AeAe/A⋆A⋆, the pressure ratio required to establish the subsonic-supersonic isentropic flow sketched above must be a very specific value. The pressure ratio is a boundary condition applied to the flow in the laboratory, it is provided by a high-pressure air reservoir at the inlet or a vacuum source at the exit.
We're trying to illustrate the application of various types of CFD techniques which is why we're intentionally choosing a flow problem with an already known analytic solution.
CFD SETUP: Here, we'll set up three echelons of equations as follows:
1. The governing flow equations will be couched in terms of partial differential equations suitable for the time marching solution of quasi-one-dimensional flow.
2. The finite-difference expressions pertaining to MacCormack's technique as applied to this problem will be set up.
3. Other details for numerical solution(such as the calculation of the time step and the treatment of boundary conditions) will be formulated.
We know that the governing equations can be formulated in two forms:
1. Conservative form 2. Non-conservative form
When we are applying the quasi-one-dimensional assumptions to these equations, then we get the equations suitable for the present problem.
The governing equation for quasi-one-dimensional are as follows:
1. NON-CONSERVATIVE FORM:
Continuity Equation:
A∂p∂t+V∂(ρA)∂x=-ρA∂V∂xA∂p∂t+V∂(ρA)∂x=−ρA∂V∂x
Momentum Equation:
ρ∂V∂t+ρV∂V∂x=-∂p∂xρ∂V∂t+ρV∂V∂x=−∂p∂x
Energy Equation:
ρ⋅Cv∂T∂t+ρVcvρ⋅Cv∂T∂t+ρVcv∂T∂x∂T∂x= -p∂V∂x-pV∂(InA)∂x−p∂V∂x−pV∂(InA)∂x
We're now going to proceed by setting up our numerical solution for the above set of equations. Now, we're aware of these terms being dimensional variables. This has an added engineering advantage as it lets us have a feel for the magnitude of the actual magnitudes of these physical quantities as the solution progresses eventually. However, for nozzle flows, the flow field variables are being referenced to their reservoir values where the dimensional flow variables are divided by their respective reservoir values.
The non-dimensional variables P/P0P0, ρρ/ρ0ρ0and T/T0T0 vary between 0 and 1, which is an "aesthetic" advantage when presenting the results as fluid dynamicists dealing with nozzle flows so frequently use these non-dimensional terms, we will follow suit here.
A number of CFD practitioners prefer to always deal with non-dimensional variables, whereas others prefer dimensional variables, as far as the numerics are concerned, there isn't much real difference and the choice is really a matter of our personal preference.
Now, here reservoir temperature and density are denoted by T0T0and P0P0, respectively, we define the non-dimensional temperature and density, respectively as:
T' = T/T0T0 ρ′ρ'= ρρ0ρρ0
Moreover, the length of the nozzle is being denoted by 'L' and hence we define our dimensionless length as:
x′=xLx'=xL
Denoting the speed of sound in the reservoir as 'ao' we define a dimensionless velocity as:
V′=VaoV'=Vao
Also, the quantity L/a0a0has the dimension of time, and we define a dimensionless time as:
t′=tLaot'=tLao
Finally, we're considering the ratio of the local area A to the sonic throat area A⋆A⋆ and dimensional pressure divided by reservoir pressure 'P0P0'.
We define a dimensionless area and pressure respectively as:
A' = A/A⋆A⋆ P' = P/PoPo
Then returning to the governing equations and introducing these dimensionless variables, we get the non-dimensional governing equations as follows:
CONTINUITY EQUATION: ∂ρ′∂t′=-ρ′∂V′∂x′-ρ′V′∂(InA′)∂x′-V′∂ρ′∂x′∂ρ′∂t′=−ρ′∂V′∂x′−ρ′V′∂(InA′)∂x′−V′∂ρ′∂x′
MOMENTUM EQUATION: ∂V′∂t′=-V′∂V′∂x′-1γ(∂T′∂x′+T′ρ′∂ρ′∂x′)∂V′∂t′=−V′∂V′∂x′−1γ(∂T′∂x′+T′ρ′∂ρ′∂x′)
ENERGY EQUATION: ∂T′∂t′=-V∂T′∂x′-(γ-1)T′[∂V′∂x′+V∂(InA′)∂x′]∂T′∂t′=−V∂T′∂x′−(γ−1)T′[∂V′∂x′+V∂(InA′)∂x′]
DISCRETIZATION OF EQUATIONS: Finally, we've set up a particular form of equations that will be most appropriate as well as convenient for the time-marching solution of quasi-one-dimensional nozzle flow. This method consists of the simplified form of solving PDE's.
As second-order degree equations are hard to solve, by the use of this method, we get second-order accurate time marching without much effort. So, the Macormack method is an easy form of explicit finite-difference technique which is second-order-accurate in both space and time. It consists of two steps, Predictor and Corrector.
We're starting with the Predictor Step. In this step, we set up the spatial derivatives as forwarding differences. Also, to reduce the complexity of the notation, we'll drop the use of prime to denote a dimensionless variable. What follows is, all variables are non-dimensional variables, denoted earlier by the prime notation. As it is an explicit method, we find the time derivative terms from the previous state values and update them in the flow field variables. Then we proceed with the corrector step. We're then calculating the time derivatives by using these updated values by using backward difference. After that, the average values of these time derivatives are found out.
By using these average values, we have to find the corrected values of the flow field variables at time t+dt. These steps should be repeated until a steady-state is reached.
CALCULATION OF TIME-STEP(ΔtΔt):
Δt=C(Δx√T+V)Δt=C(Δx√T+V)
Then, we've to calculate ΔtΔtat all grid points ranging from i=1 to i= n while taking C as '0.5' & then choosing the minimum value among them. The resulting ΔtΔt obtained is used in the above equations to calculate flow field variables.
NOZZLE SHAPE & INITIAL CONDITIONS:
The nozzle shape, A = A(x) is specified and held fixed, independent of time, for the case illustrated in this section, we choose parabolic area distribution given by:
A=1+2.2(x-1.5)2A=1+2.2(x−1.5)2 0≤x≤30≤x≤3
Hence, for our choice of initial conditions, they're chosen close to steady-state values so as to reduce the chances of the solution going unstable.
In the present problem, we know that ρρ and T decrease and V increases as the flow expands through the nozzle. Thus, we are choosing the initial conditions that qualitatively behave in the same fashion. For simplicity, we are now assuming linear variations of the flow-field variables, as a function of x. For the present case, we assume the following values at time t = 0.
Initial Conditions at t = 0,
ρ=1-0.3146⋅xρ=1−0.3146⋅x
T=1-0.2314⋅xT=1−0.2314⋅x
V=(0.1+1.09⋅x)T12V=(0.1+1.09⋅x)T12
BOUNDARY CONDITIONS:
Boundary conditions are important because the applicability of numerical methods and the resultant quality of computations can critically be decided on how those are numerically treated and without the proper implementation of boundary conditions and their numerically proper representation, we have no hope whatsoever of obtaining a proper numerical solution to our flow problem.
Many important problems such as flow driven by moving objects, free-surface flow, flow involving air bubbles, flow accompanying phase transition, and fluid-structure interaction, are moving boundary problems. In dealing with such boundaries with movement or deformation, the traditional mesh methods such as finite difference method, finite element method, and finite volume method generally encounter difficulties in calculating accurately the geometric shape of the boundary.
SUBSONIC INFLOW BOUNDARY AT POINT 1:
Here, we're allowing one variable to float & we're choosing the velocity V1 because on a physical basis we know that the mass flow through the nozzle must be allowed to adjust to the proper steady-state and hence allowing variable with velocity 'V1' to float makes the most sense as part of this adjustment. The value of V1 changes with respect to time and the information is obtained from the flow-field solution over the internal points. We're using linear extrapolation from point 2 & 3 to calculate V1 and we take ρ(1)ρ(1)= 1 & T(1)= 1 as fixed values.
V(1)=2⋆V(2)-V(3)V(1)=2⋆V(2)−V(3)
SUPERSONIC OUTFLOW BOUNDARY AT POINT N:
For this, we must allow all flow-field variables to float. We're then choosing to use linear extrapolation based on the flow-field values at the internal point. Specifically, we have the values for the non-dimensional variables.
VNVN= 2VN-12VN−1-VN-2VN−2
ρNρN= 2ρN-12ρN−1-ρN-2ρN−2
TNTN= 2TN-12TN−1-TN-2−TN−2
Now, we have to write a program to use this data to perform simulation in order to find out the flow field variables and required data.
2. CONSERVATIVE FORM:
CONTINUITY EQUATION:
∂(ρA)∂t+∂(ρAV)∂x=0∂(ρA)∂t+∂(ρAV)∂x=0
MOMENTUM EQUATION:
∂(ρAV)dt+∂(ρAV2)∂x=-A∂p∂x∂(ρAV)dt+∂(ρAV2)∂x=−A∂p∂x
ENERGY EQUATION:
∂(ρ(e+v22)A)∂t+∂(ρ(e+V22)AV)∂x∂(ρ(e+v22)A)∂t+∂(ρ(e+V22)AV)∂x
We then convert these equations into a non-dimensional form where we get the governing equations as:
CONTINUITY EQUATION:
∂(p′A′)∂t′+∂(ρ′A′V′)∂x′=0∂(p′A′)∂t′+∂(ρ′A′V′)∂x′=0
MOMENTUM EQUATION:
∂(ρ′A′V′)∂t+∂(ρ′A′V2+(1γ)p′A′)∂x′=1γp′∂′A∂x′∂(ρ′A′V′)∂t+∂(ρ′A′V2+(1γ)p′A′)∂x′=1γp′∂′A∂x′
ENERGY EQUATION:
∂[ρ′(e′γ-1+γ2V′2)A′]∂t′+∂[ρ′(e′γ-1+γ2V′2)V′A′+p′A′V′]∂x′=0∂[ρ′(e′γ−1+γ2V'2)A′]∂t′+∂[ρ′(e′γ−1+γ2V'2)V′A′+p′A′V′]∂x′=0
Lets express the equations for the quasi-one-dimensional flow in generic form. We're going to define the solution vector U, the flux vector F, and the source term J as follows:
U1U1= ρ′A′ρ′A′
U2U2= ρ′A′V′ρ′A′V′
U3U3= ρ′(e′γ-1+γ2V′2)A′ρ′(e′γ−1+γ2V'2)A′
F1F1= ρ′A′V′ρ′A′V′
F2F2= ρ′A′V′2+1γρ′A′ρ′A′V'2+1γρ′A′
F3F3 = ρ′(e′γ-1+γ2V′2)V′A′+p′A′V′ρ′(e′γ−1+γ2V'2)V′A′+p′A′V′
J2J2= 1γp′∂A′∂x′1γp′∂A′∂x′
Now, we're substituting these terms in the above equations to get,
∂(U1)∂t′=-∂F1∂x′∂(U1)∂t′=−∂F1∂x′
∂(U2)∂t′=-∂(F2)∂x′+J2∂(U2)∂t′=−∂(F2)∂x′+J2
∂(U3)∂t′=-∂(F3)∂x′∂(U3)∂t′=−∂(F3)∂x′
Now, to obtain primitive variables(P,V,T,ρρ,etc) we must decode the elements U1,U2 & U3 as follows. From the definitions of U1, U2 and U3 given above, we have:
ρ′ρ′= U1U1 / A′A′
V′V′= U2U2/ U1U1
T′T′= e′e′ = (γ-1)(γ−1)(U3U1-γ2V′2)(U3U1−γ2V'2)
p′=ρ′T′p′=ρ′T′
In the above equations, we can notice that the flux vector elements F1, F2 & F3 are source terms couched in terms of the primitive variables. In the program when F1, F2 & F3 are expressed directly in terms of ρρ, P, V, e, instabilities develop during the course of the time-marching solution. So, we write such that the governing equations are in pure form in terms of the elements of the solution vector(In terms of U1, U2 & U3).
Finally, we're now going to do discretization. We're taking the governing equations that consist of flux terms and source terms and discretizing them by using MacCormack's Method.
PREDICTOR STEP: To start the calculation, we're using the initial conditions for U1, U2, and U3 to calculate the initial values of F1, F2, and F3 & we use these values to calculate the values of U1, U2, and U3 in the next time step.
Calculating the time derivatives by using the forward difference to the RHS terms to discretize and these time derivatives are used to calculate the predicted values of U1, U2, and U3.
CORRECTED STEP: Now using the predicted values, we can calculate the new flux terms and source terms. We're then calculating the time derivative terms by using backward difference & then we're calculating the average values of time derivatives of predictor and corrector steps.
After that, we're then calculating the corrected values of U1, U2, and U3 by using the average values of time derivatives. This is done at every node except the boundary nodes since we calculate them by using boundary conditions. Then we're also calculating the primitive variables. These steps are to be repeated until a steady state is obtained.
INITIAL CONDITIONS:
The nozzle's shape, A = A(x), is specified and held fixed, independent of time. For the case illustrated in this section, we're choosing a parabolic area distribution given by:
A=1+2.2(x-1.5)2A=1+2.2(x−1.5)2 0≤x≤30≤x≤3
Dependent variables which are being solved from the governing equations are U1, U2 & U3. We need initial conditions for these same variables at time t = 0 in order to start the finite-difference solution. The initial conditions for U1, U2, and U3 will also allow initial conditions for F1, F2 & F3 to be obtained from the above equations. Such initial conditions for F1, F2 & F3 are needed to form the x derivatives on the right-hand sides of the governing equations at the first time step.
The initial conditions for U1, U2, and U3 were synthesized by assuming the following variations for ρ′ρ′ & T′T′:
These variations are slightly more realistic than those assumed for the nozzle shape & initial conditions for non-conservative form is made with such anticipation that the stability behavior 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 given in case of non-conservative form.
The initial condition for the variation of V′V′ is synthesized by taking advantage of the fact that one of the dependent variables in our governing equations, namely U2 is physically the local mass flow; that is U2 = ρ′A′V′ρ′A′V′. Therefore, for initial conditions only, let's assume a constant mass flow through the nozzle and calculate V′V′ as
V′=U2ρ′A′=0.59ρ′A′V′=U2ρ′A′=0.59ρ′A′
The value 0.59 is chosen for U2 because it is close to the exact analytical value of the steady-state mass flow(which for this case is 0.579). Therefore initial condition for V′V′ as a function of x′x′ is obtained by substituting the p′p′ variations given by the above equations. Finally, the initial conditions for U1, U2 & U3 are obtained by substituting the above variations for p′p′.
T′T′& V′V′ into the definitions above:
U1U1 = ρ′A′ρ′A′
U2U2 = ρ′A′V′ρ′A′V′
U3U3 = ρ′(e′γ-1+γ2V′2)A′ρ′(e′γ−1+γ2V'2)A′
BOUNDARY CONDITIONS:
SUBSONIC INFLOW BOUNDARY CONDITIONS AT POINT 1:
At the subsonic inflow boundary, two properties are held fixed and one is allowed to float & at the supersonic outflow boundary all properties are allowed to float. In the present formulation, just as before, we're holding p′p′and T′T′ fixed at the inflow boundary, both equal to 1.0 and it allows V′V′to float. By holding p′p′ fixed, then the U1 at grid point i =1 is fixed independent to time, via U1 = p′A′p′A′. That is:
U1(i=1)U1(i=1) = (ρ′A′)i=1(ρ′A′)i=1 = (A′)i=1(A′)i=1 = fixed value
U2(i=1)U2(i=1) = 2U2(i=2)2U2(i=2)- U2(i=3)U2(i=3)
U3U3 = U1U1(T′γ-1+γ2V′2)(T′γ−1+γ2V'2)
SUPERSONIC OUTFLOW BOUNDARY AT POINT N:
In the case of the supersonic outflow boundary, all properties are allowed to float. The flow properties at the downstream, supersonic outflow boundary are obtained by linear extrapolation from the two adjacent internal points. If N denotes the grid point at the outflow boundary, then:
(U1)N(U1)N= 2(U1)N-1(U1)N−1 - (U1)N-2(U1)N−2
(U2)N(U2)N = 2(U2)N-1(U2)N−1 - (U2)N-2(U2)N−2
(U3)N(U3)N = 2(U3)N-1(U3)N−1 - (U3)N-2(U3)N−2
MATLAB CODE TO SOLVE THE QUASI 1D NOZZLE FLOW PROBLEM USING BOTH CONSERVATIVE AND NON-CONSERVATIVE FORMS:
1. NON-CONSERVATIVE FORM:
clear all
close all
clc
% Inputs
n = 31;
x = linspace(0,3,n);
dx = x(2)-x(1);
% Calculate Initial Profiles
rho = 1-0.3146*x;
t = 1-0.2314*x; % temperature
v = (0.1+1.09*x).*t.^0.5;
% Area
a = 1+2.2*(x-1.5).^2;
nt = 1400;
dt = 0.04;
gamma = 1.4;
th = round(n/2);
cfl = max((v+sqrt(t))*dt/dx);
% Outer Time Loop
for k = 1:nt
rho_old = rho;
v_old = v;
t_old = t;
% Predictor Method
% Continuity Equation
for j = 2:n-1
dvdx = (v(j+1)-v(j))/dx;
drhodx = (rho(j+1)-rho(j))/dx;
dtdx = (t(j+1)-t(j))/dx;
dlogadx = (log(a(j+1))-log(a(j)))/dx;
% Continuity Equation
drhodt_p(j) = -rho(j)*dvdx - rho(j)*v(j)*dlogadx-v(j)*drhodx;
% Momentum Equation
dvdt_p(j) = -v(j)*dvdx-(1/gamma)*(dtdx+(t(j)/rho(j))*drhodx);
% Energy Equation
dtdt_p(j) = -v(j)*dtdx-(gamma-1)*t(j)*(dvdx+v(j)*dlogadx);
% Solution Update
rho(j) = rho(j) + drhodt_p(j)*dt;
v(j) = v(j) + dvdt_p(j)*dt;
t(j) = t(j) + dtdt_p(j)*dt;
end
% Corrector Method
for j = 2:n-1
dvdx = (v(j)-v(j-1))/dx;
drhodx = (rho(j)-rho(j-1))/dx;
dtdx = (t(j)-t(j-1))/dx;
dlogadx = (log(a(j))-log(a(j-1)))/dx;
% Continuity Equation
drhodt_c(j) = -rho(j)*dvdx - rho(j)*v(j)*dlogadx-v(j)*drhodx;
% Momentum Equation
dvdt_c(j) = -v(j)*dvdx-(1/gamma)*(dtdx+(t(j)/rho(j))*drhodx);
% Energy Equation
dtdt_c(j) = -v(j)*dtdx-(gamma-1)*t(j)*(dvdx+v(j)*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 i = 2:n-1
rho(i) = rho_old(i) + drhodt(i)*dt;
v(i) = v_old(i) + dvdt(i)*dt;
t(i) = t_old(i) + dtdt(i)*dt;
end
% Apply the boundary conditions
% Inlet Boundary Condition
v(1) = 2*v(2) - v(3);
% Outlet Boundary Condition
v(n) = 2*v(n-1) - v(n-2);
rho(n) = 2*rho(n-1) - rho(n-2);
t(n) = 2*t(n-1) - t(n-2);
%Pressure
p = rho.*t;
%Mach no
M = v./sqrt(t);
%Mass flow rate
m = rho.*a.*v;
figure(3)
hold on
if k == 1
plot(x,m);
end
if k == 50
plot(x,m);
end
if k == 250
plot(x,m);
end
if k == 500
plot(x,m);
end
if k == 750
plot(x,m);
end
if k == 1400
plot(x,m);
end
rho_throat(k) = rho(th);
v_throat(k) = v(th);
t_throat(k) = t(th);
p_throat(k) = p(th);
M_throat(k) = M(th);
end
xlabel('Nozzle length')
ylabel('Mass flow rate')
legend('1dt','50dt','250dt','500dt','750dt','1200dt')
title('Variation of mass flow rate at different time steps')
figure(1)
subplot(5,1,1)
plot(x,v)
xlabel('Nozzle length')
ylabel('Velocity')
title('Steady state primitive variable profiles')
subplot(5,1,2)
plot(x,rho)
xlabel('Nozzle length')
ylabel('Density')
subplot(5,1,3)
plot(x,t)
xlabel('Nozzle length')
ylabel('Temperature')
subplot(5,1,4)
plot(x,p)
xlabel('Nozzle length')
ylabel('Pressure')
subplot(5,1,5)
plot(x,M)
xlabel('Nozzle length')
ylabel('Mach Number')
figure(2)
subplot(5,1,1)
plot(v_throat)
xlabel('Time step')
ylabel('Velocity')
title('Transient variation of primitive variables at throat')
subplot(5,1,2)
plot(rho_throat)
xlabel('Time step')
ylabel('Density')
subplot(5,1,3)
plot(t_throat)
xlabel('Time step')
ylabel('Temperature')
subplot(5,1,4)
plot(p_throat)
xlabel('Time step')
ylabel('Pressure')
subplot(5,1,5)
plot(M_throat)
xlabel('Time step')
ylabel('Mach Number')
figure(3)
plot(x,m);
title('Steady state mass flow rate profile')
xlabel('Nozzle length')
ylabel('Mass flow rate')
2. CONSERVATIVE-FORM:
clear all
close all
clc
% Inputs
n = 31;
x = linspace(0,3,n);
dx = x(2)-x(1);
% % Calculate Initial Profiles
% rho = 1-0.3146*x;
% t = 1-0.2314*x; % temperature
% v = (0.1+1.09*x).*t.^0.5;
for i = 1:n
if x(i) <= 0.5
rho(i) = 1;
t(i) = 1;
end
if x(i)> 0.5 && x(i)<= 1.5
rho(i) = 1 - 0.366*(x(i)-0.5);
t(i) = 1 - 0.167*(x(i)-0.5);
end
if x(i)> 1.5
rho(i) = 0.634 - 0.3879*(x(i)-1.5);
t(i) = 0.833 - 0.3507*(x(i)-1.5);
end
end
nt = 1400;
dt = 0.0267;
gamma = 1.4;
th = round(n/2);
% Area
a = 1+2.2*(x-1.5).^2;
v = 0.59./(rho.*a);
U1 = rho.*a;
U2 = rho.*a.*v;
U3 = U1.*(t/(gamma-1) + (gamma/2)*v.^2);
cfl = max((v+sqrt(t))*dt/dx);
% Outer Time Loop
for k = 1:nt
U1_old = U1;
U2_old = U2;
U3_old = U3;
F1 = U2;
F2 = U2.^2./U1 + ((gamma-1)/gamma)*(U3 - (gamma/2)*(U2.^2./U1));
F3 = gamma*U2.*U3./U1 - 0.5*gamma*(gamma-1)*(U2.^3./U1.^2);
% Predictor Method
% Continuity Equation
for j = 2:n-1
% Continuity Equation
dU1dt_p(j) = -(F1(j+1)-F1(j))/dx;
J2 = (1/gamma)*rho(j)*t(j)*(a(j+1)-a(j))/dx;
% Momentum Equation
dU2dt_p(j) = -(F2(j+1)-F2(j))/dx + J2;
% Energy Equation
dU3dt_p(j) = -(F3(j+1)-F3(j))/dx;
% Solution Update
U1(j) = U1(j) + dU1dt_p(j)*dt;
U2(j) = U2(j) + dU2dt_p(j)*dt;
U3(j) = U3(j) + dU3dt_p(j)*dt;
end
% Predictor update
rho = U1./a;
t = (gamma-1)*(U3./U1 - (gamma/2)*(U2./U1).^2);
F1 = U2;
F2 = U2.^2./U1 + ((gamma-1)/gamma)*(U3 - (gamma/2)*(U2.^2./U1));
F3 = gamma*U2.*U3./U1 - 0.5*gamma*(gamma-1)*(U2.^3./U1.^2);
% Corrector Method
for j = 2:n-1
% Continuity Equation
dU1dt_c(j) = -(F1(j)-F1(j-1))/dx;
J2 = (1/gamma)*rho(j)*t(j)*(a(j)-a(j-1))/dx;
% Momentum Equation
dU2dt_c(j) = -(F2(j)-F2(j-1))/dx + J2;
% Energy Equation
dU3dt_c(j) = -(F3(j)-F3(j-1))/dx;
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 i = 2:n-1
U1(i) = U1_old(i) + dU1dt(i)*dt;
U2(i) = U2_old(i) + dU2dt(i)*dt;
U3(i) = U3_old(i) + dU3dt(i)*dt;
end
%Final update
rho = U1./a;
t = (gamma-1)*(U3./U1 - (gamma/2)*(U2./U1).^2);
v = U2./U1;
% Apply the boundary conditions
% Inlet Boundary Condition
U1(1) = a(1);
U2(1) = 2*U2(2) - U2(3);
U3(1) = U1(1)*(t(1)/(gamma-1) + 0.5*gamma*v(1)^2);
% Outlet Boundary Condition
U1(n) = 2*U1(n-1) - U1(n-2);
U2(n) = 2*U2(n-1) - U2(n-2);
U3(n) = 2*U3(n-1) - U3(n-2);
%Pressure
p = rho.*t;
%Mach no
M = v./sqrt(t);
%Mass flow rate
m = rho.*a.*v;
figure(3)
hold on
if k == 1
plot(x,m);
end
if k == 50
plot(x,m);
end
if k == 250
plot(x,m);
end
if k == 500
plot(x,m);
end
if k == 750
plot(x,m);
end
if k == 1400
plot(x,m);
end
rho_throat(k) = rho(th);
v_throat(k) = v(th);
t_throat(k) = t(th);
p_throat(k) = p(th);
M_throat(k) = M(th);
end
xlabel('Nozzle length')
ylabel('Mass flow rate')
legend('1dt','50dt','250dt','500dt','750dt','1200dt')
title('Variation of mass flow rate at different time steps')
figure(1)
subplot(5,1,1)
plot(x,v)
xlabel('Nozzle length')
ylabel('Velocity')
title('Steady state primitive variable profiles')
subplot(5,1,2)
plot(x,rho)
xlabel('Nozzle length')
ylabel('Density')
subplot(5,1,3)
plot(x,t)
xlabel('Nozzle length')
ylabel('Temperature')
subplot(5,1,4)
plot(x,p)
xlabel('Nozzle length')
ylabel('Pressure')
subplot(5,1,5)
plot(x,M)
xlabel('Nozzle length')
ylabel('Mach Number')
figure(2)
subplot(5,1,1)
plot(v_throat)
xlabel('Time step')
ylabel('Velocity')
title('Transient variation of primitive variables at throat')
subplot(5,1,2)
plot(rho_throat)
xlabel('Time step')
ylabel('Density')
subplot(5,1,3)
plot(t_throat)
xlabel('Time step')
ylabel('Temperature')
subplot(5,1,4)
plot(p_throat)
xlabel('Time step')
ylabel('Pressure')
subplot(5,1,5)
plot(M_throat)
xlabel('Time step')
ylabel('Mach Number')
figure(3)
plot(x,m);
title('Steady state mass flow rate profile')
xlabel('Nozzle length')
ylabel('Mass flow rate')
OBSERVATIONS :
STEADY-STATE DISTRIBUTION OF PRIMITIVE VARIABLES INSIDE THE NOZZLE:
1. STEADY-STATE VALUES OF FLOW FIELD VARIABLES IN NON-CONSERVATIVE FORM:
2. STEADY-STATE VALUES OF FLOW FIELD VARIABLES IN CONSERVATIVE FORM:
The above plots are almost the same when compared to the analytical solution referred from Anderson Textbook. We can observe from the above graphs that while the Velocity & Mach Number increase from inlet to outlet, it comes at an expense of a drop in pressure, density & temperature along the same path of the inlet to outlet. From Bernoulli's equation we know that in the case of horizontal flow, an increase in velocity must be accompanied by a decrease in pressure because of conservation of energy. The energy required to increase the fluid's velocity comes at an expense of static pressure energy.
Here, the most important parameter to understand is Mach Number. Mach Number increases with an increase in value of x.
The Mach number is a dimensionless quantity. Mathematically, the Mach number equation can be written as: M = V / C , such that V is the speed of the object and C is the speed of sound.
VARIATION OF MASS FLOW RATE DISTRIBUTION INSIDE THE NOZZLE AT DIFFERENT TIME STEPS DURING THE TIME MARCHING PROCESS:
TIME-WISE VARIATION OF MASS FLOW ALONG WITH THE AXIS IS PLOTTED HERE USING THE NON-CONSERVATIVE FORM:
Here, the non-dimensional mass flow is plotted as a function of non-dimensional distance through the nozzle at different time steps in the graph during the time-marching procedure. The first two curves represent the variation of profile at initial conditions. We can see that after 250 iterations, the flow achieves an almost constant mass flow profile.
TIME-WISE VARIATION OF MASS FLOW ALONG WITH THE AXIS IS PLOTTED HERE USING THE CONSERVATIVE FORM:
1. Conservative & Non-Conservative form has different initial conditions, due to which we get a different mass flow rate curve when compared to the non-conservative form analysis.
2. The curves are closer to each other and also to their analytical value which is 0.579 which means that more mass is preserved in the case of the conservative form than the non-conservative form.
FIGURE SHOWING THE TIME-WISE VARIATION OF VARIOUS PRIMITIVE VARIABLES AT THE THROAT USING THE NON-CONSERVATIVE METHOD:
FIGURE SHOWING THE TIME-WISE VARIATION OF VARIOUS PRIMITIVE VARIABLES AT THE THROAT USING THE CONSERVATIVE METHOD:
From the above plots, we can see that the time-wise variation of non-dimensional flow field variables at the throat in both non-conservative form and conservative form respectively attaining the steady-state. We can also observe that the oscillations are reduced to a steady-state eventually.
Here, we can see that after almost 200 timesteps, the non-conservative form attains a steady-state while the conservative form takes around 300 timesteps for it to attain a steady state as the oscillations continue for a bit longer in its case. Therefore we can conclude that the non-conservative form attains steady-state faster than the conservative form.
CONCLUSIONS DRAWN FROM OBSERVATIONS ABOVE:
1. The conservation form does yield a better mass flow distribution. It is better at preserving the mass throughout the flow field mainly because the mass flow itself is one of the dependent variables in the equations as the mass flow rate is the primary result derived from the equations of conservation-form. The dependent variables in the non-conservation form of the equations are primitive variables, and the mass flow is obtained only as a secondary result in this case.
2. Even though the final solution for both forms is similar, the non-conservative forms oscillate more when compared to the conservative form.
3. Non-conservation form leads to smaller residuals and the amount by which the residuals decay is often used as an index of the quality of the numerical algorithm. Hence, we know that the non-conservation form does a better job.
4. No clear superiority of either form is observed in terms of accuracy of the results although the non-conservation form seems to produce slightly more accurate results for the flow field variables/ primitive variables and the conservation form seems to produce slightly more accurate results for the flux variables.
5. Judging by the amount of effort needed to achieve a solution we can conclude that the conservation form requires more effort most of which is 5ue to the fact that we need to decode the primitive variables from the flux variables and such decoding is not necessary when we're solving for the problem using the non-conservation form.
6. Both the forms: conservative and non-conservative are grid-independent.
Thus, MacCormack Method provides quite stable and close to exact solutions hence, it is widely implemented to solve various flow problems.
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...
FINAL GD&T PROJECT: BUTTERFLY VALVE WITH GD&T IN SIEMENS NX CAD
OBJECTIVE: The primary objective of this project is to design and model individual components of a butterfly valve using the provided drawings while applying Geometric Dimensioning and Tolerancing (GD&T) principles to each component within the Siemens NX CAD environment. Upon successfully creating the individual…
13 May 2024 10:55 AM IST
WIRING HARNESS FLATTENING & DRAWING WORKBENCH
OBJECTIVE: Take the harness assembly from the previously completed challenge and flatten it. Position this flattened view on the drawing sheet. It’s important to make sure that bundles with protective coverings are visually distinct in the drawing view. This step is part of our ongoing process to create a drawing…
13 May 2024 09:30 AM IST
FINAL PROJECT TWO: BACKDOOR WIRING HARNESS USING CATIA V5
OBJECTIVE: This project aims to demonstrate the practical application of wiring harness routing and design principles on a car's backdoor/tailgate using CATIA V5 software. The main objective is to showcase the implementation of industry best practices and packaging rules studied throughout the course by creating a properly…
15 Apr 2024 07:58 AM IST
FINAL PROJECT ONE: V16 ENGINE WIRING HARNESS ROUTING, PACKAGING, FLATTENING AND DRAWING
OBJECTIVE STATEMENT: The primary objective of this assignment is to design and route a comprehensive wiring harness for a given engine using CATIA V5 software. The design process will encompass applying industry-standard packaging rules, best practices, and guidelines acquired through the coursework. Particular emphasis…
08 Mar 2024 06:46 AM IST
Related Courses
Skill-Lync offers industry relevant advanced engineering courses for engineering students by partnering with industry experts.
© 2024 Skill-Lync Inc. All Rights Reserved.