Tip:
Highlight text to annotate it
X
We have seen how the CFD method works in principle for the case of finding out the laminar flow
steady, fully developed flow through a duct of rectangular cross-section. As we have seen
the method, let us now try to work it out and get to know it better how the CFD solution
works. So, we are going to take again the case of a rectangular cross section. We will
solve the corresponding Gaussian equation for it with actual numbers to see how the
solution works using the CFD approach.
So, what we are looking at is steady fully developed laminar flow through a straight
duct of rectangular cross section. So, we take the geometry to be like this. The equation
that we are interested in is dou square w by dou x square plus dou square w by dou y
square, equal to d p by d z times 1 by mu. When we want to work out, we have to put some
numbers to this. So, we will put this whole thing to be equal to 1. This is the equation
that we want to solve subject to this domain and we will put some numbers to this very
soon and subject to the condition that w is equal to 0 on all walls.
So, the first step of the CFD solution is that we have a governing equation and we have
the boundary condition. We have a flow domain with actual dimensions of this flow domain
yet to be specified and we should put here this as being equal to minus 1 because mu
here is the dynamic viscosity which is constant, and pressure gradient is negative in the positive
flow direction. So, if you want w to be positive then we should put a negative sign so, we
will just put for the nominal case of this value on the right hand side being minus 1.
The second step of the CFD solution is to discretize the flow domain and identify the
grid points at which we want to evaluate the solution. This is a hand calculation. So,
in order to make this possible, we will take small number of grids. We will divide this
into four parts, four equal parts and this into three equal parts, just to show the unequal
thing. We say that, this is where we measure; we want to find out the velocities. As per
the boundary condition that w is equal to 0 on all walls; these are the wall and points;
this is where the velocity is 0. So, in this we would like to know the velocities
at these six interior points subject to the condition that w is equal to 0 at all the
other grid points, and that the resulting w distribution within this at these grid points
satisfies this relation approximately. So, we have to put some number here. So, we will
put i here in the x equal to constant planes to be 0 1 2 3 4 and j equal to 0 1 2 3.
So, in this problem the unknowns can now be listed. These are the value of w at these
six points and we can identify them with the i index and the j index. So, these are w i
comma j and since both i and j are in single digits we can write this as w i j. Specifically,
we are interested in this point which becomes w 1 1, this is w 2 1, this is w 3 1, this
is w 1 2, this is w 2 2 and this w 3 2. So, there are six unknowns and these are the values
that we want to find subject to the w distribution satisfying this condition, and also this equation
where the right hand side is given as minus 1.So, the pressure gradient divided by viscosity
for this fully developed flow within this rectangular duct is this.
Now, here we are also taking the case of delta x being constant and let us put this as one
unit. At this stage, we do not have to specify what that unit is, may be meters or millimeters
and so on. It depends on the actual problem. For the sake of the numerical example, we
will take delta x to be 1 and delta y to be 1. We will do later on an assignment where
these numbers are different. So, under these conditions, now everything is specified here
for us to launch our third part of the CFD problem.
So, we have taken the equations, we have discretized the points and we have located the unknowns
here. So, this is step one, step two, and third step is to write down the approximate
form of this equation which is going to be valid at each point, and as we have said,
we will approximate dou square w by dou x square, at the location i comma j. As w i
plus 1 j, minus 2 w i j, plus w i minus 1 j, divided by delta x square, and dou square
w by dou y square at i comma j is approximated as w i j plus 1, minus 2 w i j plus w i j
minus 1 by delta y square. So, this equation is going to be discretized. Using these approximations
as; this plus this is equal to minus 1 at a point i j.
So, we have finally the discretized equation. We can also take a note of the fact, that
delta x is equal to 1 and delta y is equal to 1. So, we can take a note of this. So,
the denominator is equal to 1 here and is also equal to 1 here. So, we can directly
substitute this. As the discretized equation at point i comma j is given by this, plus
this equal to minus 1, i comma j plus, w i minus 1 j plus, w i j plus 1 minus, 2 w i
j plus, w i j minus 1 is equal to minus 1. So, this formula we apply at each of the points
at which we want to have the solution. So, for example we take this recurrence, this
formula here. We apply at this point, this point, this point and that’ll give us six
equations that we can get from this. So, for this point we can say that point 1 comma 1
is this one where it is 1 in this direction, 1 in this direction. So, this means that i
is equal to 1, j is equal to 1. So, we substitute this. Here we get, w 2 1 minus, 2 w 1 1 plus,
w 0 1 1 plus, w 1 2 minus, 2 w 1 1 plus, w 1 0 is equal to minus 1.
So, these points with zeros are the boundary values. For example, w 0 1 is this point and
w 1 0 is this point and w 0 1 is this point. and we know from the boundary condition that
at all these point; w is equal to 0. So, where we know the boundary value, here we can substitute
this so, this is equal to 0 and this is equal to 0, and we have minus 2 w 1 minus 4 w 1
1. So, we can write this as minus, 4 w 1 1 plus w 2 1, plus w 1 2 equal to minus 1. We
take this thing to the other side, so, we can get w 1 1. We take it here and then we
divide by minus 4. So that’ll be, 1 plus w 2 1 plus w 1 2 divided by 4. This is the
formula for w 1 1. We will do the same thing for the next point
which is this one. So, that is point 2 1. This implies that i is equal to 2 and j equal
to 1. So, we substitute these values again in this equation. We should be getting w 3
1 minus 2 w 2 1, plus w 1 1, plus w 2 2, minus 2 w 2 1, plus w 2 3. This is 2 3 and this
is 2 1 and this whole thing is equal to minus 1. Here, there are no boundary points, so,
we can write this as minus 4 w 2 1 is equal to, we can take all the other things onto
the right hand side minus w 3 1, minus w 1 1,minus w 2 3. So, this is 2 0, this is j
is equal to 1, so this is 0 here and 2 0 is this one and this is 0 here. So there is a
wall point. So this is like this and so from this we can get w 2 1 is equal to 1 plus w
3 1, plus w 1 1, plus w 2 3 divided by 4. We can do the same for this point; that is
3 1 this implies i is equal to 3 and j is equal to 1.So, let us see, if we can do without
mistake. So, w 4 3 almost so, w 4 1, minus 2 w 3 1, plus w 2 1, plus w 3 2, minus 2 w
3 1 and this is w 3 0.
So again w 3 0 is the point at this value, and w 4 1 is this point. This is 0 and this
value is 0.So, this is 0 here and this is 0. If these values are different, we can substitute
those values, and then take them to the right hand side. So, if the boundary condition is
known, then that goes onto the right hand side, it would not appear on the left hand
side as the equation. So, we can rewrite this equation, as minus 4 3 1, plus w 2 1 plus
w 3 2 equal to minus 1. So, this gives us, w 3 1 as 1 plus w 2 1. So, this is the equation
for w 3 1 and we can similarly do for point 1 2. Here, this implies i is equal to 1 and
j is equal to 2. So, this is w 2 2 minus, 2 w 1 2, plus w 0 2, plus w 1 3, minus 2 w
1 2, plus w 1 1 is equal to minus 1. We are looking at this, so we have 0 2 is 0 here
and we should be getting 1 3 here. So, we should be getting 2 wall values here, so this
is 0 and this is 0. So, 1 3 is 0 and 0 2 is 0, so both are wall values.
So, from this we get the expression for w 1 2 1 plus w 2 2 by 4. So, this is the expression
for w 1 2. So, we have to do these two expressions and if we want, we can look at a pattern that
is developing. So, the value at this particular thing is given as the right hand side here,
plus the neighboring points this plus this, plus this, plus this, these two are 0 divided
by 4, we can go through that. So, for example, 1 plus w 2 1, plus w 3 2 by 4, so we can do
it like this or we can do it the hard way. Let us continue the hard way, so this is 2
2 1 is equal to 2 and j is equal to 2. Substitution of this expression into this, w 3 2, minus
2 w 2 2,plus w 1 2,this w 2 3, minus 2 w 2 2, plus w 2 1 is equal to minus 1. So, we
are looking at this point here, so we have w 2 3 is a wall value. So, this is equal to
0, so it will be 1 plus this, plus this, plus this, divided by 4.
So, from this we can write w 2 2 is equal to 1,plus w 3 2,plus w 1 2,and the last point
is w point 3 2 and here i is equal to 3,and j is equal to 2. So, substituting this, here
we get w 4 2 minus 2 w 3 2, plus w 2 2, plus w 3 3, minus 2 w 3 2,plus w 3 1 is equal to
minus 1. Here, we are considering this point, so this is an average of these four. So, there
should be one boundary point here and one boundary point here, so, that is 4 2 and 3
3. So, this is 0 and 4 2 is 0. So, we can write from here, w 3 2 is equal to 1, plus
w 2 2, plus w 3 1 divided by 4. So, by approximating at each of these discrete points at which
we want to get the solution, by approximating that partial differentiation equation with
finite difference formulas, and coming up with a template here after substitution of
delta x and delta y and the constants like this, we can take this template and apply
it to each point, at which we want to get a solution, and from that template we will
be able to derive the corresponding algebraic equations which give us the values. We are
expressing the value at w 1 1 here in terms of these points, these expressions like this.
So, we find that there are six equations, a, b, c, d, e, and f which can be solved.
These are algebraic equations and these are such that we cannot solve any of these equations
individually. So, that is, we cannot from this equation, find w 1 1 and this equation
or from this equation w 2 1, like that. So, these are simultaneous algebraic equations
and these are also equations with constant coefficients. All of them have a coefficient
of one, or zero or something depending on what we have here, and we can see that we
have a set of linear equations with a linear algebric equations with constant coefficients,
which we get as a result of discretizing the governing partial differential equation at
each point, at which the value or the variable is to be evaluated, so we have converted.
If you put all these into a w equal to b form, then the matrix a here will have the coefficients
of all the variables and we can write them down, and we can then invert the matrix. So,
that is one way of doing it. But, in this example, typically in computational fluid
dynamics, we do not directly invert the matrix for a solution of a w equal to b. It is because,
we normally deal with very large matrices, so inversion of those matrices is very expensive
computationally, so we use some iterative methods. So, we illustrate the application
of the Gauss-Seidel iterative method. It is a fairly basic iterative method for the solution
of a of i equal to b. At this stage, we take it as granted that this method will work for
this set of equations. We will take a minute to explain, what this, how this iterative
method works and then we can actually look at the solution.
So, we consider a set of three equations involving three variables, a 1 1 x 1, plus a 1 2 x 2,
plus a 1 3 x 3 is equal to b 1. a 2 1 x 1, plus a 2 2 x 2, plus a 3 3 x 3 is equal to
b 2. And a 3 1 x 1, plus a 3 2 x 2, plus a 3 3 x 3 is equal to b 3. So, these are a system
of three linear algebraic equations involving variables x 1, x 2, x 3 with constant coefficients
a 1 1 to a 3 given here, and also b 1, b 2, b 3 are given here, and we want to find out
from this number x 1 x 2 x 3.
There are many ways of doing this and what we have in our example here is the set of
six equations involving six variables and we can write down all the six in the form
of a 6 by 6 equation like this, and in the Gauss-Seidel method we solve this iteratively
by expressing taking in the first equation. We express x 1 in terms of x 2 and x 3, and
we take this equation x 2 in terms of x 3 and x 1, and then x 3 in terms of x 1 and
x 2 using the third equation and we start with some guess value. We proceed with the
iterative update of x 1 x 2 x 3 from those given values until we reach a solution, which
is not changing for everywhere. So, this iterative procedure is very easy
to implement, it will work only under certain conditions. We will look at the conditions,
for convergence of this particular method. For the time being we just describe here,
how the method works and how to implement it and to get a solution. From this, so from
this we can write this as x 1 equal to b 1, minus a 1 2 x 2, minus a 1 3 x 3 divided by
a 1 1. x 2 is b 2, minus a 2 1 x 1,minus,this is a 2 3 a 2 3 x 3 divided by a 2 2, and this
is written as x 3 equal to b 3, minus a 3 1 x 1,minus a 3 2 x 2 divided by a 3 3.
So, we have written these expressions like this and we start with an initial guess which can in the case of a convergent method.
In cases where the Gauss-Seidel iterative method works, it does not matter what the
initial guess is, the final iterative solution will work in all cases. So, we can start with
an initial guess 0 0 0 for x 1 x 2 x 3 and we can get from this x 1 x 2 x 3 at the end
of the first iteration. By successively substituting in this, and then from this we get again x
2 x 3 in the second iteration by again substituting. And then we go on x 1, x 2, x 3 for the kth
iteration, until we reach successive values of x 1, x 2, x 3 which are not changing.
So, the way that we evaluate this is that, x k plus 1. At the end of k iterations, we
have the value of x 1, x 2, x 3. So, the new value of x k plus 1, is evaluated as x k plus
1is b 1 minus a 1 2 x 2 k, and a 1 3 x 3 k and x 2 is evaluated as b 2 minus a 1 2 x
1. We already have the k plus oneth value. So, this is k plus 1 and this is k here. x
3 is evaluated as b 3 minus a 3 1 x 1 k plus 1 and k plus 1 here. So, in that way, given
the kth values of x 1,x 2, x 3, we can evaluate x 1 k plus oneth value of x 1 x 2 x 3 like
this and in this way we generate from a starting guess. We can generate values and we can get
the solutions
So, this is the method that we will use for the evaluation of this and we can therefore
write these equations. Already, these are put in this Gauss-Seidel form, so we can say
that w 1 1 k plus 1 is w 2 k. Here, this is k plus 1 is 1 plus k and we already have k
plus 1, so we can put k plus 1. Here this is k. Here 3 1 k plus 1 is w 2 1. We have
k plus 1, w 3 2 we do not have yet. So, this has k w 2 k plus 1 is 1 plus w 2 2 k w 1 k
plus 1 w 2 2 k plus 1. 3 2 we do not yet have, so this is at kth value. This is what we already
have this. So this is k plus 1 and this is also k plus 1. And finally the value of w
3 2 is obtained as, k plus 1, k plus 1. So, we have here the formula for evaluating w
1 1, w 2 1, w 3 1, w 1 2, w 2 2 and w 3 2, from known values of the same variables at
the previous iteration. So, we start with an initial guess, w i j equal to 0. And then,
we generate solution at the end of the first iteration, the second iteration, third iteration,
like that. Then, as we keep on going, we will see how we can get, how we will be getting
converge solutions.
So, this iterative procedure is quite common and is best illustrated when we tabulate the
results. So, we will try to tabulate the result, so as to get an understanding of how the solution
evolves in this particular case with using this particular iterative method.
In evaluating this point, we made a small mistake here. When j is equal to 1, so this
is 2 i j plus 1, so this is w 2 2, and this 2 1, and this 2 0. So, this we put correctly,
this we put correctly here we put as 2 3, so that has to be corrected here and also
here. So, this is the formula. Now, we will try to do the solution of this using the Gauss-Seidel
method in Excel because it is very easy to see and we can also see how the solution evolves.
So, we will take a look at this we will try to program the solution of these using the
Gauss-Seidel method in Excel. So we would like to get back it
Sir there is a mistake! There is a mistake.
You want me to do now. Are they going to see that or...
Only this sir, only problem. So, just we will go into this. Here we have
the formulas for the six points w 1 1, 2 1, 3 1, 1 2, 2 2, 3 2 which are written as per
what we have derived. For example w 1 1 is 1 plus w 2 1 plus w 1 2 divided by 4. Like
this, and here we are associating the six variables and this is the solution, at the
initial guess. So, we are giving this as zero values, for each of this and the value here
w 1 1 is coded as per this 1 plus w 2 1, is this one here. So, that is, b 2 plus w 1 2,
is d 2 divided by 4. And if you, consider for example 1 2 is given as 1 plus w 2 2,
so that is this value plus w 1 1 and since we take the latest value it should be this.
So it is, 1 plus e two plus a 3, divided by 4. So, that is how we have programmed. So
this is, for w 1 1. This is w 2 1, which has w 3 1, plus w 1 1, plus w 2 2.
So, these are coming up like this w 3 1. We have already seen that is w 2 1 for which
we have the latest value and w 3 2 for which we have the previous value here. w 1 2, w
2 2 has w 3 2 which is yet to be evaluated. So, that is F 2 w 1 2 which is already evaluated,
so it becomes D 3. This one here and w 2 1 which is already evaluated which is this B
3 divided by 4.And finally w 3 2, here is w 2 2 which is evaluated, so we put it as
this value that is E 3 and w 3 1 which is already evaluated. So, that is E 3 by 4. So,
we have put the formulas here and then let us just erase all this.
So, this is what the value is at the end of the first iteration and we can take this down.
We can just click and the solution will evolve. So, we can see this is the iteration number,
so we can put this as iteration number. This is equal to the first iteration that is the
initial guess. This is the second iteration, so we can write this as 1 plus 1 plus, so
second iteration. We have these values, third iteration, these are changing and we can see
that as we go more and more iterations it is this particular w 1 1 has gone from 0.25
to 0.406, 0.497, and 0.543. But, it is beginning to stabilize by the end of the sixteenth iteration.
It has become stable up to the fifth decimal point. Here, similarly is these also have
become stabilized, and we have got the final solution which is accurate up to the fifth
decimal place as given by this. So, this is how, we have generated the solution for the
six unknowns using the Gauss-Seidel method. Using an iterative scale, if you were to change
the initial guess, for example, if you make this as 1, so if you change everything to
1 the solution is changing but the final solution will not change. The final converge solution
will not change, I can even make this as 10 here, the solution is not changing. The converged
solution will not change if I change my initial guess. So, that is one characteristic feature
of a convergent iterative method. What will change of course, is the values as they go
through the iterative process and the number of iterations which are needed to get to the
converge solution may change but the final solution itself will not necessarily change.
So, in this way we can get a solution to the overall scale. So, this in a sense, the CFD
approach. What we have illustrated here is for the very simple case with just six points,
and six points will not give us a very good solution. We have to do it for much larger
number of points so that we have to do much more work, but it will be done. We are going
to show the final solution that we get, when we increase the number of grid points and
from what we have here, we have four divisions in this direction and three divisions in this
direction. We can make it bigger like this; there are
more number of grid points in this direction. More number of grid points and obviously even
these values become unknown. Now, the unknowns are as many as this and we can make it even
more. By putting further, subdividing in this way and in this way, so the more the number
of grid points, the more accurate will be the solution. We will see how the solution
changes when we put this. When we subdivide it further and further at some point we will
notice that with further subdivision the solution will not change, and that is the point at
which we have a grid independent solution. That is the ultimate solution to the governing
equation that we have started out with.
We have seen that doing out by hand is very tedious, now we come to the computer. And
then once we have a computer, we can easily program this and get a solution for large
number of grids. Here, we have a table which summarizes the solution for different grid
sizes when you have a 6 by 6, 7 by 7, 10 by 10, and 50 by 50,100 by 100, like this. We
can see that as the grid size increases the matrix size also increases. For example, if
you have 41 points by 41points then we have matrix size of 1600 by 1600. We have close
to 1600 variables, in this, and in each case, we can use the Gauss-Seidel method and then
get a converged solution. Once we have the converged solution, we can find out the average
velocity by taking an average of all the values. And that average values, what is put here
in this column and we can see that for our 6 by 6 grid size. We have an average velocity
of 0.7409, this is the average speed, and as the grid size increases to 10 by 10 it
has come down to 0.6722. For 25 by 25, it has come down to 6070 and as the grid size
increases, it seems to be stabilizing. For example at 51 by 51, it is 5843,at 71 by 71,
it is 0.5765,so the changes becoming less and less, and 94 by 94 is 0.5742, so it is
about 57 is the velocity that we are getting as the average velocity. This is what we expect
as average velocity and the same thing is plotted here. It is the mean velocity that
is computed versus the grid size. We can see that for small grid sizes, there is a rapid
change here, but eventually it is coming down to about 0.57 and what is plotted in these
three figures are the contours of the velocity. Within the cross section, we can see that
at for a grid size of 15 by 15, we have high velocity at the center and otherwise we have
almost like an elliptic type of things showing different velocity levels and at the four
corners we have low velocity, as you increase the size here it is stabilizing to a particular
pattern of the velocity. So, here you can see rough edges. Whereas
here with the increased number of grid points, we have smooth contour lines which is again
an indication that reaching a stable solution stable solution in the sense, a converged
solution which is known as a grid independent solution, So, this illustrates, the power
of computational fluid dynamics, that we can now get a velocity variation within the cross
section by having large number of grid points spread throughout. We can use fairly simple
methods which are very easy to program in the computer and which are therefore easy
to solve and then get the velocity distribution for this fairly complicated shape under steady
conditions. So, the methods that we have employed in this
are I would say elementary and which requires only elementary knowledge of mathematics.
It does not require much more than this. So, there lies the appeal of computational fluid
mix that the underlying principles are something that can be understood fairly easily by an
engineer, with a basic course in fluid mechanics and mathematics.
So, what we have done in this lecture is to put in practice the CFD method of our solution.
We have taken the partial differential equation which governs the flow and we have arrived
at a template of an approximate form of the partial differential equation which is applicable
for any grid point within the flow domain. We have used this template to generate as
many number of algebraic equations as there are unknown variables within representing
the velocity at the interior points of that particular flow domain.
This gave us a matrix, a w equal to b. The size of the matrix depends on how many grid
points that we take, and we have specifically used one method which is useful for solving
large number of these equations which is the Gauss-Seidel method. If you apply that, then
go through the iterations preferably and almost invariably on a computer, and then we get
a converged solution. The converged solution can be integrated or added up to compute the
average velocity. We can find out how the average velocity changes, with the increase
in grid size, and ultimately get a grid independent average velocity, at which point we can plot
the contours, like what we have shown here. This represents the variation of velocity
within the interior, so this how a CFD method would work. But we have deliberately taken
a very simple method and the governing equation is very simple. The geometry of the flow domain
is very simple, the discretizing method is also very simple and the solution method is
also very simple. CFD can do much more than what is shown in
this simple example and this is what we are going to do in the next several classes. We
will introduce complexity, in each of these and then come up with a method which can ultimately
solve general three dimensional flow which is unsteady, which may be turbulent, which
may be reacting and which may be three dimensional. So, in an arbitrarily complex geometry, we
will try to come up with a solution method which is robust efficient and which can also
consider the realistic geometries that are present in industrial equipment and realistic
flow conditions that are present in industrial equipment and using these techniques one should
be able to look at a number of scenarios, like what is the optimal distribution of a
mixing equipment within a reactor, so as to get a highest mixing efficiency at the lowest
power consumption in a chemical reaction engineering application. Also, if you are looking at a
boiler furnace, what kind of firing practice must be adopted for a given coal in order
to minimize the formation of nox gases. So, if you are looking at those kinds of realistic
scenarios, one has to be able to solve three dimensional unsteady reacting turbulent flows.
The principle is that we have enunciated here in this will be extended to deal with those
kinds of complexities in the coming lectures. So, in the next lecture we are going to start
with an outline of the overall course having now known what CFD is, and then we will go
into the each of the modules which represents complications on each of these things.