Tip:
Highlight text to annotate it
X
In the last lecture, we have seen one particular way of dealing with the flow domain which
does not fit into the usual coordinate system that we that we normally adopt, for example,
the Cartesian or cylindrical and all that. And if we have a flow domain which is a combination
of these, for example, a pipe inside a rectangular duct then we cannot make use of this. And
in the last lecture, we have seen an approach whereby we come up with a new coordinate lines,
which distort and wrap around the bodies of interest so as to define the flow domain in
a structured coordinate frame: psi, eta, zeta instead of x, y, z.
And we solve the equations, we transform the equations from the physical plane into the
computational plane and we map the physical plane into the computational plane, and then
we do all the computations that are the discretization, solution and all that. In the computational
plane and then map, the solution back into the physical plane. So, this body fitted approach
is one way of dealing with complicated geometry and it has proved very successful and it retains
the characteristics of a structured mesh approach to the solution of the governing equations.
Now in this lecture, we will see in an alternative view point, an alternative approach. It is
based on the finite volume method; essentially what we are saying is that we want to solve
our governing partial differential equation to get the solution of the flow variable at
several points within the flow domain that is a basic CFD approach. And these partial
differential equations represent the conservation in terms of the rate of increase of that particular
property in a control volume as being balanced by three causes: one is the convective flux,
the other is the diffusive flux and the third one is the source or the sync term or the
source and the sync term in the case. For example, of the turbulent kinetic energy,
where we saw that there is a production term associated with a large eddies large eddies
and then there is a sync term associated with the with the very small eddies and together
they determine, the overall level of the particular quantity, which is being conserved; which
is the turbulent kinetic energy, in this case at that particular grid point in that particular
control volume.
So, we can rewrite our standard conservation equation, which we have written many times
as in this way. This is the convective flux is equal to the diffusive flux plus the source
term. So, this particular equation can be rewritten in this way by bringing the fluxes
together of rho u j phi minus gamma phi d phi by d x j equal to S phi.
So, this is the convective flux and this is the diffusive flux of the particular quantity phi. And in turbulent flow, we have seen that this particular term
will dominate and this diffusibility here is not that of molecular diffusibility, but
of turbulent diffusibility. But generally, this is the form here and which can be written
realizing that this is a vector divergence operator here. You can write this as plus
del dot, for example, F equal to S phi, where F is the total flux. We can put F phi here,
where F phi is F phi convective plus F phi diffusive and the convective flux is given
by rho u j phi and the diffusive flux is given by rho phi like this.
So, this is athis is a governing equation which is valid at every point. In the finite
volume method, we put the equation in this in this form. And we take a control volume
and then we integrate this over the control volume dou by dou t of dou phi dv plus integral
of del dot F phi dv equal to S phi dv. So, for every control volume that is for we construct
a volume around each grid point at which we want to evaluate the variable.
So, we can say that this is integrated over a control volume and a control volume is obviously
a closed surface. This is the control volume at the different sides and making up together
is the total volume contained in this and it has a closed surface and taking advantage
of it. We can using Gausses law, we can convert this volume integral into a surface integral.
And we can write this as dou by dou t of rho phi plus integral over the control surface
the closed control surface of n ,where n is the normal vector associated with the surface
enveloping the control volume equal to s phi d v, where this is the source term here.
So, this is a rewriting of the governing equation in what can be called as a conservation form
and the conservation form refers to the fluxes that are appearing here, the convective flux
and the diffusive flux. The conservation form refers to these terms that are appearing here
being interpreted as the convective flux and the diffusive flux and therefore associating
these with the with the fluxes that are passing through the faces of the control volume, which
make up which make up the overall volume of that.
So, giving the meaning the interpretation for these two terms as convective and diffusive
fluxes and evaluating these things over the surface enveloping the control volume is what
we call as the conservation form. And for example, the alternative view is that you
have rho u j phi and even if you take rho to be constant, we know that u is not constant
and we can write this equivalently as rho u j dou phi by dou x j plus phi dou by dou
x j of rho u j. And so these two are mathematically equivalent and similarly, we can write this
as gamma phi dou square phi by dou x j square plus dou phi by dou x j times dou gamma phi
by dou x j. So we can split it up into that and when we
do that then we lose the interpretation of this term as the flux, essentially what we
are saying is that the quantity phi for which we are writing the conservation equation;
it will the value of this within this control volume will change. That is what this term
is? This term will be non zero either, because we have a source term which is spread throughout
the control volume or because we have some flux, which is coming into the control volume
by diffusive action by the fact that the value of the point of the variable in this control
volume is less or more than the value of the control volume in the neighboring cells or
by the fact there is a flow, because of the flow and because the flow brings in all the
fluid properties. And phi is one of the fluid properties like enthalpy or temperature or
concentration. So, the value of the phi within this control
volume may change, because it is being brought in and taken out along with the flow or because
it is been diffused by gradients that exist between at this particular point, at this
particular cell across the… So, these are the mechanism by which the phi can change
and that is what is being represented in this interpretation.
Now, when we write in this interpretation, it becomes easy for us to apply this statement
of the conservation, of the conservation law here to an arbitrary control volume. Not necessarily,
something that has four faces in two dimensions or six faces in two dimensions. We can even
take a triangular control volume and then we can say that this this and we can for example,
if you say that this is A, B, C here. We can say that dou by dou t of rho phi. Let us say
that this is the ith cell with three triangular phase rho rho phi i plus…
We can associate with this the control surface three phases AB, BC, CA and over each of this
we can write. So, we can say that n dot F phi ds over AB plus n dot F phi times ds over
side BC plus n dot F phi ds over side CA, that is the integral of all the fluxes that
are coming through the phases and there are three phases here is equal to S phi times
dv, we have a dv here. So, in this particular case, the volume we are looking at the two-dimensional
thing we can make it into three-dimensional by considering unit length in the in the other
direction. So, we need to know what the volume of this
element is, so that is we need to the area of this triangle multiplied by the unit distance.
So, this is given by the cell geometry and in evaluating this term this is also given by the cell geometry and orientation. And
similarly these things and here we have again cell geometry and what we need in order to
apply this convert this equation into a mathematical equation is obviously, we have d rho phi by
dt times volume. So, for a given cell we know this delta volume
and for a given cell with control phases we know each other areas of each of the side
AB, BC and CA. So, when we talk about the area of the side. For example, BC you are
talking about the area of this surface with this side as BC and this delta z direction as unity as the length here.
So, that is how we get surface and an area and this surface has an outward novel vector
n, which is what is coming here.
So, in order to apply this equation, we need to evaluate the fluxes. So, the the flux term
F p has to be evaluated. So that means that flux coming from the convective and flux coming
from diffusive contributions will have to be evaluated at each phase. And that convective
flux and diffusive flux has a particular mathematical form rho u j phi and then the diffusive flux
and this is where we need information of the values of the phi in the neighboring cells.
So, using those values and using the definition of these fluxes here. We can come up with
the evaluation of these terms, in terms of the neighboring cell values. And once you
put this here, we will have an overall equation for this particular cell. So, that cell will
then be in the form of… At that point, we can write dou phi by dou t as rho phi i n
plus 1 minus rho phi i n minus 1 by delta t times delta volume is this and the fluxes
will have contribution associated velocity and the value of the phi here and that is
value of the phi at… For example, at this point that is the cell centroid here, the
plane the surface centroid and that has to be evaluated from in by some means.
So, we can say that this plus net convective plus diffusive flux over all surfaces enveloping
surfaces is equal to the s phi at i times delta volume. So, this is the discretized equation in which
one can see that phi i n is coming phi n plus 1 is coming and obviously in the diffusive
fluxes, we have to evaluate the gradient dou phi by dou x and that means that phi i plus
1 minus phi i minus 1 by two delta x is one possible form. So that is where we are going
to link the value of the current cell value to the neighboring cell values.
So, in this way, we can take a governing equation; we can divide our flow domain into cells,
small cells and over each cell we can discretize the conservation equation written in conservation
form. That is the variation within the cell as a result of volumetric forces and the fluxes
here, and it is in this form, we can apply this to a particular cell and then convert
this differential equation into an algebraic equation. Now, there are when we talk about
the corresponding momentum equation here. We have to make a distinction between sources
here that are truly volumetric and those are in a way coming through the surface.
So, when if when you treat the volumetric source as a volume and then as a volumetric
source term and any source term which is coming, which is acting on the control phases as appropriately
as the phase related source term then we have a strong conservation form. I would like to
say that in the momentum equation, the source term here is for example, minus dou p by dou
x is something is one of the source term. And so, we can say that minus dou p by dou
x times the volume delta v is the momentum source that is coming from the pressure, but
we know that pressure is the surface force. So, this should be actually evaluated as the
as something like equivalent to a pressure acting on the surface of the control volume
and then that would make it, what is known as strong conservation form.
Again we have… If we have radiation as a source, radiation source term is typically
associated with radioactive flux. We can so this is also a flux term, this should be coming
through the surfaces not as a volumetric source term, but if you were to evaluate this as
a volumetric source term then it is not in the strong conservation form.
So in the evaluation of source terms, we should try to make a distinction between those sources
that are coming through acting through the phases like the pressure and the radioactive
flux and those that are truly volumetric source terms like a gravitational force term is a
volumetric or is a source term which is acting throughout the material which is there in
the in the volume. If you have, for example, for the energy equation; if you have a heat
generation term, which is paid throughout the control volume then it has to be treated
as a volumetric source term, but energy flux coming through the phases in the form of radioactive
flux must be treated in in this. So, we have to treat each term appropriately
as either acting through the phases of the control volume or having or being spread out
through the throughout the domain and therefore constituting a volumetric source term. So,
once we do that then we can claim to have a strong conservation form and the conservation
equation that is being put out here, that is being solved exactly at each control volume
and that is one of the strong point of of the finite volume method.
The conservation that is incorporated in this equation is being enforced in each control
volume by interpreting this as a conservation equation in this conservation form. And so,
if the flux terms are evaluated properly, then there is no possibility of generating
spurious mass sources or flux sources or sources of this and that is one of the advantages
of this. And when you have discontinuities in the within the flow domain as may be arising
from shocks in those kind of things. Evaluating in this conservation form is supposed
to give superior results than putting it in the general form like this. So that, enforcement
of the conservation in each each discrete cell is a characteristic future of this and
written in this way this can be applied to an arbitrary control volume with defined volume
and a defined phases which envelope and completely close the cell. So, this is the basic idea,
we will we will try to see this in an action through a simple example, before we go onto
the formal evaluation of these fluxes and then look at some more complexities.
So, the case that we are looking at the simple example, that we are looking at is flow through
fully developed flow fully developed laminar steady flow through a triangular duct, through a duct of triangular cross section.
Why did we take this particular example, because this is one can see that at once straight
away; it is difficult to get an analytical solution. If it is a circular pipe, we can
do and if it is a rectangular pipe may be with more difficulty we can do, but if it
is an arbitrarily shaped triangular duct then to get the velocity field is not trivial thing.
Analytically and also to apply this to… For example, to fit a structured grid for
this is again not an easy task and definitely we will have large distortions of the cells,
around this corner points. So, if you were applying the body fitted grid
approach for this then one can expect more difficulty with the numerical solution, and
that is also one reason, why this particular example is a suitable example to illustrate
the benefits and the advantages of the finite volume method. So, in when we take this particular
cross section, we obviously need to have the governing equation and the governing equation
is the one dimensional momentum equation through the z direction.
So, and it is since the flow is steady and fully developed. It takes a form of dou square
w, where w is the velocity through this duct. So, this is the duct and velocity w is in
this direction, and x and y are within the plane of this plus dou square w by dou y square
times mu is equal to minus 1 by rho dou p by dou z. And dou p by dou z is constant for
the case of fully developed flow and it is a given constant from the fact that you can
specify the boundary condition at z equal to zero and z equal to l ,because the you
have pressure gradient as constant, we can evaluate this. So, this is a given constant.
So for this particular case we can look ,we can see ,how we can apply the finite volume
method for evaluation of this, and we can see how we are going to tackle this and evaluate
each of the fluxes through this example.
So, here we have the case of fully developed laminar flow in a triangular duct of irregular
cross-section. When we say irregular cross-section, the three sides are not the same length; so,
it is not an equilateral triangle to make it slightly more complicated and the governing
equation for this is obviously dou square w by dou x square plus dou square w by dou
y square equal to 1 by nu dou t by dou z. And we are considering constant properties
and since it is fully developed and steady this right hand side is a constant, we are
calling it as C. And the boundary conditions are, because we have three walls here the
velocity w is 0 on all the walls. So, that is the boundary conditions, so the problem
statement is very straight forward, but analytical solution is not available for an arbitrary
triangle. We have analytical solutions for equilateral
triangles which are given in stand books, but we can use this readily using the finite
volume method.
So, as a precursor, we write the governing equation in conservation form. So, what we
have here as this term, we recognize that this is the diffusion term and therefore it
is a divergence term. So, we put this in the divergence term. So, that part associated
with the left hand side here is expressed as a divergence, that is del dot gradient
of w is equal to c. And this divergence form when integrated will give us; we can convert
this into the area integral, surface integral. So, we divide the domain this triangular domain
into a combination of triangles and rectangles. And that is the advantage of finite volume
method, it is not necessary that the domain has only rectangles or only triangles; it
can be a combination of all these things. And that makes it an unstructured mesh and
we are cheating a bit here. We are taking, for example, we are dividing this into cells
like this, a combination of rectangular cells and triangular cells. And we are saying that
we want to have the value of the w at these points. These are obviously the centre centroids
of these rectangles and we are for the sake of simplicity, we are saying that we want
to have the velocity at this point, at these points here. And we know that these points
are on the wall, therefore the velocity there is 0.
As we make it more and more, as we introduce more and more number of points that approximation
will become less and we can also become more sophisticated and put this point to be at
the centroid of this triangle, but for the time being we will say that this is where
we want to have the velocities. So, the problem reduces to finding the velocities at these
4 points plus 6 points each of which has a control volume or an area which is unequal
here. And so we… The governing equation is now
integrated of the control volume and is converted into a surface integral using Gauss’s Divergence
theorem. So, that is del dot gradient of w dV integration over the control volume, closed
volume is equal to n dot gradient of w this is that diffusive flux of momentum in this
particular case, over the control surface. And that is equal to the constant right hand
side times, the integrated over the same volume. So, if we can apply this to each cell in this
way, so that is integral over the control surface is equal to this n is is an outward
normal vector for this case of two-dimensional case. It has two components i direction, S
x in the i direction, S y in the j direction. And the gradient of velocity also has two
components dou w by dou x in the i direction and dou w by dou y in the j direction. So
the dot product of these two times the constant C times, the volume of that control volume
is the discretized form of this equation. So, each cell gives an algebraic equation
linking the cell value that is w i, we can see obviously, the w i coming from the gradients
here with those of the neighboring cells.
So, we are looking now at the practicalities of the of the of the domain. And here we are
looking at 20 nodes; we are breaking up into the left hand side of this triangle into 4
equal parts here, and then 4 equal parts on this thing. And based on that we can do decomposition
of this into a combination of triangular and rectangular tiles and again the right hand
side of this is also made into 4 divisions in the horizontal and 4 divisions in the vertical.
So, we have a total of 20 points: 1,2,3,4,5,6,7,8,9,10,11,12 all the way like this, out of which some points:
1,9,15 and all those things lie on the boundary. So, out of the 20 points 8 are boundary points
with 0 velocities. The velocity is already known and we need to find the velocity at
the other 12 points. So that, we can evaluate and having found the velocity at each point,
we can multiply for example, the velocity at 16 by the area of the cell, this cell to
get the flow rate through this. We take this velocity times this area will
give you the flow rate through this, velocity at 17 time this area will give you the flow
rate through this and then we can compute the overall flow rate in that way. And so
that, we can get the flow rate for a given pressure gradient, so that solve the problem
is post. And in each case we need to find out, the delta the sides of these phases,
so we have delta x. So, we are taking for whatever reason in this
example, this length to be 0.01375 and this length to be 0.02625. So, this is broken up
into 4. So that, you have a delta x here and this is again broken up into 4 and this whole
height is taken as 0.01452. So, that is broken up into 4. So, using these things for any
cell here rectangular or triangular; we know the sides and that is the geometric information
that we have and once we know the side for each of these. We can find out the length
of each side, which constitute the surface and the area of each cell which in three dimensional
would be the volume. So, the geometrical information is obtained
from the discretization from the breaking up of the entire flow domain into tiles, rectangular
tiles and rectangular and rectangular tiles. So, we need to find out velocity at 2, 3,
4, 5, 6, 7 and then again 10 and 11, 12, 13 and then 16 and 17. So, there are 12 unknowns
and over each of these points, we have a control volume which is defined here. And for each
cell, we apply this discretized form in order to derive an equation.
So, that is shown for cell 2 here, so that, this is cell 2 with which is a rectangular
phase and here we have w 2. And immediately to the right at the centre of the adjacent
cell, we have w 3. And the bottom wall, this bottom phase here is a wall, and left side
we have this as part of the wall here; and it has a w 1 of 0 and above this above this
again we have 0 here. So, when you look at the 4 neighboring points here. On this side
you have 3, w 3 which is to be determined. On the left hand side, we have w 1 at a certain
distance which can be evaluated. And that is 0, above this you have w 9 at a certain
distance which is 0 and below this, you have a wall. So, you do not have a point here.
So, we can say that the gradient dotted with area of each side is equal to the C by mu
times, the total volume of each cell. So, that is what we have volume of… So, we can
say that all the cell to surface integral of this over the cell to surface is equal
to C by mu times volume of cell 2. So, this particular thing we are putting it as A, B,
I, H and so the integral of this over A, B, I, H. So, that is gradient of w times dA over
AB plus gradient of w times dotted with dA over BI plus gradient times dA over IH and
gradient times this over HA is equal to this. And we take advantage we take note of the
fact that these are oriented in x and y directions the outward norm vectors are located in the
x and y directions. And this one on this phase, the outward norm vector for this is located
in the positive x direction and here it is located in the negative x direction. Therefore,
gradient of w dotted with dA on AB is in the negative x direction. So, we have minus dou
w by dou y dotted with area of this cell. So, this is delta x and in the z direction
we have delta z. And over BI this is the gradient at this point, so that is dou w by dou x component
here, and the area of this cell is delta y times delta z which is in the other direction.
And here it is the gradient in the y direction at this point times the area of this cell
which is delta x, delta z, and here we have a minus, because outward normal vector is
in the negative x direction. So, we have dou dou w by dou x which is the
gradient to the x direction times the area, which is delta y, delta z and this equal to
this. So, now in this equation, we know delta x and delta z for this particular cell. And
we need to evaluate the gradients at each of these points, now when come to this particular
point minus dou w by dou y at AB. So, we take the centre of here we know that this value
is 0 from the boundary condition and this value is w 2. So, we can say that is w 2 minus
0 divided by this distance which is delta y by 2.
So, that is the estimate of the gradient that is coming here, when you come to this point
here, the gradient dw by dx can be evaluated as w 3 minus w 2 divided by this distance.
So, w 3 minus w 2 divided by delta x is the total thing, and here again we have this point
and this point. So, w 9 minus w 2 divided by this height which is delta y w 9 is 0.
So, we can make use of this and here it is w 2 minus w 1 divided by delta x. So, we are
substituting, we are making estimates for the fluxes at each of the phases. In this
case, we are looking at the diffusive flux and we can substitute like this. And where
ever we know the value 0’s here, we can substitute and then we can rewrite this and
then we can finally get an expression an algebraic equation involving w 2 and w 3, because this
side and this side, this side the velocities are known.
So, we have an algebraic equation, by applying the conservation equation for this particular
cell and we can do the same thing for each of the other cells. For example, if you consider
this cell here, you would have the gradient to be evaluated here, here and here and here.
So, this gradient will be in terms of w 10 minus w 3, this gradient will w 4 minus w
2 3, this gradient will be w 3 minus w 2 and this gradient is w 3 minus 0 divided by this
half distance. So, in that way we can evaluate for each of these.
And then, we can come up with with set of 12 algebraic equations for the 12 phases,
for the 12 control volumes or tiles. And there in expressed as linear algebraic equations
involving the grid information which is coming in form of delta x, delta y, delta z that
is here like this and so you have 12 simultaneous linear algebraic equations.
And you can put them in the matrix form and we can see that you have this central diagonal
and some adjacent diagonal, there are some things that are coming here, but there are
also some other things that are coming here.
So, it is not necessary that, because of the of the way that the numbering of the cells
is made. It is no longer possible to have the kind of structured diagonals that we have
in a case of structured grid, because the neighbor if you look at 11 for this point
here is 10 and 12 in this direction is ok, but here is 16 and 4. It is not i j plus 1and
j minus 1. So, that kind of nomenclature is not possible.
So, that is why we do not exactly have a structured matrix, but it is a still a sparse matrix
and still we have a linear algebraic equation of the form Aw equal to b. And this we can
solve using Gauss-Seidel method provided that this satisfies the Scarborough criterion,
diagonal dominance condition. And it does satisfy the diagonal dominance condition,
even in the case of this unstructured grid.
And then you can get a solution at all these intermediate points from which you can draw
the contours. So, after solving this using Gauss-Seidel method we can get a velocity
field on this.
If you want more accuracy, we can divide this into more number of points, and here it is
divided into something like 72 points here. And we have each of them is smaller and you
get more points and then we can get an equation.
And here, this is an evaluation with each horizontal division into 8 by 4, more number
of cells 32, 64 like this. And the mean velocity, that is computed here and the Reynolds number
that is computed here and the size of the matrix. As you go to larger and larger matrixes
larger larger sizes, we can see that the computed velocity is becoming constant, and the Reynolds
symbols is also becoming constant. So, we are looking here at a solution for the mean
velocity for a given pressure gradient which which seems to be ok. And which is obtained
using a combination of rectangular and triangle cells, not using a structured mesh, but using
the control volume approach.
In which, we take the discretized that governing equation conservation form and then converted
into evaluation of surface fluxes and volume source terms. So, what we have selected here
is a very simple example, because we have seen the case where the flow is a fully developed.
And that means that some of the significant components of the fluxes do not appear. If
you look at our governing equation here, we only have the diffusive flux coming on the
left hand side and here is the source term. So, we do not have the convective fluxes and
that that makes grows simplification or trivialization of the problem, as we have seen here that
is appropriate for an example, but in in reality we have to consider the convective fluxes.
The fluxes that flux of a particular component coming associated with the flow, and that
usually means that that is a difficult more difficult thing to do, because the value of
that depends on the upstream value it does not depends on the local values. A diffusive
flux depends only on local values, and a convective flux depends on the velocity direction. And
that provides a bit of challenge, when we deal with that in the context of finite volume
method. So, what we will do next is to take the general
case the take the general case of conservation equation, which not only has the diffusion
and source term, as we have considered in these example, but also the temporal and the
convective terms into that. And we can see how we can put the whole thing together in
the form of a solution procedure, which will be applicable to the general case. The that
is it It is important to able to see, how we can put together the whole process of evaluating
the fluxes and the areas and the volumes. And once we have a good understanding, we
can see that the overall method is a very flexible method. It leaves us with the choice
of defining the cells or the tiles which make up the control volume. So, it gives, it offers
us a great deal of flexibility in terms of refinement of grid and refinement of a local
zone for improved accuracy. And those kinds of things can be more readily done in a finite
volume method, there in the case of structured method, but we will also see that the evaluation
of the convective fluxes puts makes it more difficult to incorporate higher order schemes
for the first order derivatives. And that is one of the difficulties associated
with a finite volume method. So, there are deftly advantages and there are also disadvantages
with the finite volume method. Just as we have them with the structured mesh in the
body fitted coordinate system, but comparing the two, I think one would say that finite
volume method is more user friendly, more easily adaptable method, more adaptable than
the body fitted grid approach. But, body fitted grid approach is more powerful in terms of
a systematic exploitation and the full scale implementation of a high degree schemes, compact
schemes, high compact higher order schemes is not so trivial in the case of finite volume
method. So, if you are looking at highly accurate
solutions probably the structured grid approach is better. But, if you are looking at user
friendliness and in terms of being able to tackle complex geometry without too much of
difficulty, then I think the finite volume method is preferable. We still have to do
lot of work, before we can say that we know how to use these methods for the general case
of a three-dimensional flow. So, that will be the subject of the next couple of lectures.