So, good morning in yesterday’s class we
looked at the method of separation of variables it is
an important analysis technique for analytically solving partial differential equations. The
idea is to split the p d into two o d’s introducing a constant called lambda square
and the other inherent assumption is solving the o d is
much easier comparing to solving the p d. Then, we
got the general solution then we had this boundary conditions, using the boundary conditions
we got the co we evaluated the constants a, b, c, d. When we were trying to evaluate the
constants several new things entered to picture, like there is a lambda and the solution is
satisfied and lambda is equal to 1 pi by l 2 pi by l and so on. These lambdas’s are
called Eigen values of the problem. And now and then finally,
you got the solution in terms of sin n pi x by
l and sin x n pi y by l, where sin x n pi in yesterday’s call was basically because
the temperature was exponentially decaying in
the y direction. As far as the sin n pi x by l is concerned
it is basically, a repeating of periodic boundary condition. In fact it is, because of this
sin n by n pi x by l that you are able to use this
orthogonal invoke the orthogonality property. That is why this periodic boundary condition
this 0 0 boundary condition is very important, thing is towards the end of the class we also
saw that if you have a problem where more than one boundary condition is 0. Then you
have to use the super position principle have 1
0 0 plus 0 0 1 0 0 you split the problem into several
sub problem and use the method of separation of variables we individually find the solution,
and then get a composite solution which is the additive sum of all the solutions.
But, now you can see that even for a two dimensional problem two dimensional steady steady
state steady conduction in a slab with constant properties a getting an analytical solution
was quiet involved. And then finally, your finally,
you have to evaluate a series it is not that you
have got an expression like e to the power of minus m I so e to the power. So, you have
to evaluate a series, and then the accuracy of
a solution depends on the number of terms you
retain in the series. So, it is not very very convenient though you we we will claim that
it is an analytical solution and so on. You will soon
realize that if you evaluate the heat flux using the
series, you require many many terms. It is not so easy to get the heat flux with, because
again the temperature has to be differentiated and
then the k d t by dx term will contain so many
terms in the series. And for certain other types of boundary condition
for example, you are varying one you have q which is a function of x, and then you have
a convective boundary condition on one side you have two other temperatures on the other
two sides and so on. And then it becomes extremely formidable to solve these problems.
And the moment you say thermal conductivity is not a constant and so on it is very difficult
to solve these problem. Therefore, numerical techniques were developed, and the growth
of the numerical techniques was also accelerated by the development of computers. And then
this availability of programming languages likes
FORTRAN and so on, which enabled the engineers to tackle this problem. Of course, now a
days you have got sweets like mat lab which can solve many of these problems, and you
have got libraries which solve o d’s and p d’s
and Fortran and all that and then there are commercial software like comsolve is very
powerful based on finite element, and fluent which based on finite volume which is based
on finite volume which can solve all these problems with considerable ease.
So, the engine behind all these software is some numerical technique. So, these numerical
techniques broadly the numerical techniques which are used in conduction in heat transfer
basically are either the finite difference method. As far as heat transfer and fluid flow are
concerned the first two are very very popular. The
third is not so popular, but still there are few there are some practitioners would like
to use the finite element method for both for both
fluid flow and heat transfer. Finite volume most of
the commercial software our fluent works on the finite volume method. So, essentially
what is the basic difference between finite difference
finite volume and finite element method? The finite difference is different from the other
two because the finite difference method is based
on the Taylor’s series approximation of a function, from that you get the first derivative
second derivative and then you find out what is the value of the nodal value of temperature
or velocity or whatever variable under question,
as a function of it is neighboring nodes. And
then once you set up a system finally, you set up a system of linear equations which
is solved by matrix inversion or gauss seidel iteration
and precede that is how a finite difference work.
A finite volume method is basically the control volume technique it is also called a control
volume method. Where you a taking control volume and find out what are the fluxes which
are entering, what are the fluxes which are leaving, any other heat generation which is
in the medium or for any heat storage or heat depletion
that is enthalpy increase or enthalpy decrease in the medium. So, it is basically
like a you are looking at the law of conservation of
energy as far as conduction problem is concerned and writing down all the terms and finding
out an expression for the value of the temperature at a particular node. So, there is no Taylor’s
series expansion in a finite volume method. And in a finite in the finite difference method
we finally, evaluate finally, we evaluate the
derivative at a point, but what is this point we have to be very careful it has a point
of area and volume. If the point does not have any
area and any volume then you cannot define stress
or something if for example, stress cannot be defined at a point and so on. So, finally,
we have to have that in mind. In the limit delta
x delta y tending to 0 that means delta x cannot
become 0 then you are in trouble. So, finite volume method is quiet natural, because it
is based on finite control volume. You do not
reduce it to inimitably small size and make it a
point. So, naturally the finite volume method can handle certain type of things.
So, the finite element method is actually a cousin of the finite volume method where
also the governing equation is integrated once. When
the governing equation is integrated once conditions like Nyman condition like dT by
dx is a constant of for example, dT by dx is equal
to 0, is very frequently encountered in conduction problems. If you have n surfaces in a
problem one of one or two the surfaces may be insulated. DT by dx is equal to 0 is very
easily handled by a finite element method, because the governing equation is integrated
once, and then you start finding out the equation
for a particular element then you find the equations for all the elements assemble all
these elements you get a system of linear equation,
beyond a certain state it becomes again like a finite volume method. So, you integrate
the governing equation once.
But, in finite difference method what you do is you take the governing equation and
start discretizing, i plus 1 j i minus 1 j we do
not do that in a finite element method, but anyway
that I 1 i minus 1 j i j everything will come in the end. Instead of using the i j notation
in finite volume and finite in finite volume
method particularly we use east west north south
here we use i minus 1 j i plus 1 j i j i j plus 1 i j minus 1.
Now in this lectures, we will just we will look at the basic principles underlying the
finite difference as well as the finite volume method.
The finite element method is beyond the scope of this course. So, finite difference what is the basis, the
basis is the Taylor’s series expansion. So, the basis
the Taylor’s series expansion now, let T be let us say that you want to evaluate T
T x i plus T is the temperature. I am instead of using
y I want to use T, because T is frequently used in
conduction, x is the special coordinate so T is a function of x it could be a function
of x y z whatever. Now for simplicity we say T is a
function of x i using Taylor’s series we want to
expand it as T is a function of x I, plus dT by dx at x i into delta x i plus plus higher order
terms. The higher order terms will become negligibly small when delta x is delta x i
is properly chosen, because it is going to be
raised to the highest powers of delta x. so, if delta x
is point naught one, this will be point naught one, this will be point naught one squared.
This will be point naught one cube and remember
at the bottom also it is two factorial three factorial everything is increasing. So, higher
order terms can be neglected even though d square T by dx squared maybe significant and
higher order derivatives may be significant. Now this equation can be rearranged for dT
by dx. So, rearranging equation 1. So, this delta x
i is small compared to x i small expansion, that is very very important it is x i plus
epsilon if x is 100 delta x should not be 75. If x is 100
delta x can be 1.5 0.1 like that it is called as small
expansion. Rearranging equation one we have yeah plus so the key point here is you are
anyway neglecting higher order term, if you say that dT by dx is equal to approximately
this, then this will be the order of the error associated
with the trunckage. So, the order of the error is directly proportional to delta x I, since
it is proportional to delta x i to the power of one it is
called a first order accurate scheme. Therefore, now if you write like this dt by
da is not 0 it is order of that is the leading order, do
not ask me sir what about other terms other terms are there, but other terms are much
smaller compared to this fellow. So, the maximum culprit
maximum error will be caused by delta x next term is delta x square by 2 factorial.
And anyway delta x itself is small. So, this is the
gives you the order of accuracy this is called the forward difference. It is first order
accurate, because delta x i please be careful, and we
are neglecting only upto the second order, we are
taking upto second order. We are taking up to second order but, when you get the dT by
dx we are dividing by delta x therefore, d square
T by dx square gets multiplied by 2 factorial it
gets so, one delta x i gets cancelled are you getting the point so it first order accurate.
So, if you want to use this first order accurate scheme for solving a problem, you have to
use more number of grids, because there will be
a severe what is called the this error is called the truncation error for the numerical scheme.
Now this is for one node in a problem in a two
dimensional problem. You may have dou t by dou x dou t by dou y d square T by dx d y
dou square d by dou y square, then each of this
dT by dx is represented as a finite difference, instead of a differential you are representing
it a finite difference when you assemble all this
you will get a set of simultaneous equations for T of i j. I can vary from 1 to 20 j can
vary from 1 to 20 depending on your grid. Then
when you are doing these operations the computer can store only a finite number of digits.
Repeatedly when you are doing gauss seidel method
then there will be a loss, there will be a loss of information because the computer will
truncate after a certain number of decimal that is called round off error.
So, the round of error will increase the round of error will increase if you are using more
number of grids. When you are using more number of grid delta x i will come down
therefore, the truncation error will decrease. So, there is an optimum grid size at which
the sum of these two errors is the total error
will be a minimum. So, if you plot the error what are these two,
a grid size in grid size is size of the grid delta x
increases, because number of number of computer operations are reduced, because only you
are solving on a 5 by 5. You are solving on a 50 by 50 truncation error increases. What
happened tell me is grid size this is not number of grid grid size is the size of delta
x. If the size of delta x increases truncation error
will increase, there is no doubt about it. I am not
saying number of grid points here, it is one by number of grid point it is inversely proportional to the number of grid. You can
plot the same thing as number of grid points then
you have to make it ulta round off and truncate. So, this is about forward difference it is
also possible to use a backward difference. Now what will be the backward difference,
T of next step into minus and plus will alternate plus this what is this now this should be
3 rearranging plus. So, I can just write it as the backward difference
is also first order accurate. So, it is called backward difference please note we are trying
to get expressions for dT by dx dT by dy and so on, because we are frequently encountering
partial differential equation. So, we want expressions for dou t by dou x dou t by dou
y more importantly we want expressions for dou
square d by dou x square into dou square d by dou y square, because only the mixed
derivative first derivative second derivative. First derivative in time second derivative
in space, with that all combinations will be
there unsteady heat unsteady heat conduction steady
state heat conduction op everything will be covered. Of course, I am deriving all these
only for the Cartesian coordinate, you will be
little more involved but, using the same procedure you can derive the finite difference from
approximations for the first derivative for both the
cylindrical as well as the spherical coordinates. Now approximation of the second derivative
that is very very important, because if you want
to solve the Laplace equation. The d square T by dx square and d square d square T by
d y square make an appearance. So, now we will
have to look at the approximation of the second derivative. So, please remember that even
the backward difference is also first order accurate
approximation. Now T of can I write like this itself or T
of i plus 1 can I write that i plus 1 or x of we will
keep as x of x of i plus 1, equal to this is 6 is based on forward differences. T of
so I am writing the Taylor series expansion for t
of x i plus 1. And I am also writing the Taylor’s series expansion of t of x i minus 1. No in
fact, now there are several ways of doing what
what do we want to do first, instead of taking the second derivative what I want to do first
is, with the same thing several things can be
done, we can call the title as central difference before going to this. Now subtracting 8 from 7 no no I want to find
first derivative. Subtracting subtracting subtracting 8 from 7 we get T x i T x i get
cancelled 2 dT by dx correct, second derivative gone third derivative is there plus. Now therefore,
minus minus therefore, dT by dx at I can
be written as T of i plus 1 minus T of i minus 1 divided by 2 delta x i but, the beauty now
is the order of accuracy is delta x i squared.
So, if then if you are evaluating the first derivative
at x i you are taking a node ahead i plus 1 minus i minus 1 divided by 2 delta x I,
but you get the second order accurate. So, if you use
a central difference only thing is the value of the
node at that the value of the variable at that node is not considered? You want to get
dT by dx at T of i but, T of i is not used in the expression
no problem. So, that is the therefore, this is
this is the central difference approximation this is the central difference approximation
for first derivative it is clearly seen that it
is second order accurate, because it is delta x delta x
square. This is the the second order accurate. If you add what happens shall we do that,
if you add adding 7 and 8 adding 7 and 8 we have t
of is equal to dT by dx is gone. And then 2 into then plus now you get get an expression
for d square dx by d square. Therefore, I think
we need a number did I put a number 9 I will call
this 10 I will call this 11. Therefore is equal to divided by 2 by delta x square, but
2 is 2 gets cancelled no no then there is a problem is
it correct, no the no no no there is a 2 here. What I
have done is T x i plus 1 plus T x i minus 1 is there, minus 2 is so the two remains
in the d square dx by d square. 2 factorial 2 factorial
will get cancelled this will be 24. So, this will be
delta x squared plus the square of square gets cancelled. It is still second order accurate,
this is central difference not forward difference
this is called the central difference. i plus 1 i
minus 1 everything is there what is the forward here, you have got forward and backward it
is its one elliptic equation so elliptic. So,
this is basically second order accurate. So, now we have got an expression for d square
x by dx square, it is called central difference. This can also be alternative derivation i j i plus 1 j i minus 1 j. so, I can take
an intermediate node i plus half j I can take here i minus
half j. So, I can get so dT by dx at half j is equal
to T of i plus 1 j minus T of i j divided by delta x, i
can get dT minus dx at i minus half j, that is T at i j minus T of i minus 1 j so this
is. Now I will say d square T by dx squared here is
dT by dT by dx at i plus 1 i plus half minus dT by
dx at i minus half divided by delta x. There is already a delta x here delta x and
delta x will give you delta x square and you will
again get T of i plus 1 j minus T of i minus 1 j T of i plus 1 j plus T of i minus 1 j
now you substitute this so you will get the same thing. But, when you do this way it is very difficult
the order of accuracy and all that, it is more
involved to get it this, but here it was straight forward is also another. So, there are two
ways of getting the same thing. Now we will go
on to the so we have got the finite difference representation of the first derivative and
the second derivative. So, therefore, we can take up a
simple partial differential equation, and try to get try to develop a finite difference
scheme, for getting the nodal value of a particular
node temperature at a particular node, in terms of its
neighboring node that is the so for everything for all the the goal of the finite difference
approach is to replace all derivatives by algebraic equation. It is a linear combination
of its neighbors that linear combination stems from
the fact that you can express everything in the
form of a Taylor’s series. So, the goal of a finite difference method
is finally, to use all the methods you have learnt in
linear algebra to solve out this system of equation. Therefore, if you have to be good
in finite difference or finite volumes you have to be
good in linear algebra. How to invert matrix whether it is worth inverting matrix more
than 100 by 100, how to accelerate the convergence of gauss siedl what is successive over relaxation
under relaxation. This is for steady for unsteady what is the time step you have to
use, number cell Fourier number, grid Fourier number, unconditionally stable, conditionally
stable explicit method, implicit method, semi implicit and explicit crank method and so
on. So, you have to be good in all this. Now now this is the first this is all the
preliminaries, mathematically how do you get the
finite difference representation. We are heat transfer engineers so our goal is to solve
the heat transfer equation. So, we will just take up
the Laplace equation and see how we can use this
to solve the Laplace equation is it clear up to this stage. d yeah yeah yeah correct
is that. Let us consider a rectangle slab same same
old, steady state 2 dimensional constant properties no heat generation. Regarding equation Del
square del square T equal to 0 no T is naught T is
naught. So, what is the equation number here, so we have to so we will give a number for
this this one there is half business we have what
is the equation number for this 13. So, we will
give a number 14 this, so d square d by dx square plus. Now, we write the finite difference expressions
for both the d square T by dx squared and d
square x d by Dy squared. For that what we will do is we will consider a node i j we
consider the 4 neighbors’ i plus 1 j and get the
second derivative with respect to x and the second
derivative with respect to y for T. any problem. So, I keeping it as delta x delta x cannot
be changing the node, but let us start with uniform
grid so number 16. Now we can add d square by let us say let
delta x be equal to delta y, let delta x be equal to
delta y therefore, d square T by dx squared means therefore, therefore, it is very simple.
So, the value of the temperature in the particular
node is basically the arithmetic average of all
the 4 nodes if which are equally which are equispaced from the node i j. so, this is
the application of the this is the application
of the finite difference method to Laplace equation.
So, it will be slightly changed if you want to apply this to the Poisson equation where
there will be a volumetric heat generation. Volumetric
heat generation can take place in a medium if there is nuclear reaction or a chemical
reaction which is taking place. For example, if you want to incorporate heat
conduction I am sorry heat generation plus
correct if q v is equal to 0 it will reduce to equation 19 equation 22. So, that q v delta x square by k v will give
an idea will give an idea of how much will be the
contribution from the heat generation term. For example, if the heat generation even if
the heat generation 1 into 10 to the power of
6 heat generation 1 into 10 to the power of 6.
Suppose if delta x is just 1 millimeter or then delta x will be 1 into 10 to the power
of minus 3, delta x square will be 1 into 10 to the
power of minus 6. So, 1 into 10 to the power 6 into 1
into the power minus 6, if you have a thermal conductivity value of minimum of 200 then
the it will it will not really matter, but if
k it is Bakelite or if it cork and then delta x is small is
larger than q v delta x x square by k v will will play it is part in the algorithm.
This can be reworked for a non uniform grid space. Non uniform grid size where delta x
i is from delta y i. so, here you will certain
T I will get multiplied by delta y by delta x and all
that. So, this essentially the formulation. So, you can just keep on formulating for other
steady equivalence for the cylindrical as well as the spherical coordinate.
In the next class first 10 minutes we will just take a 100 0 0 problem apply it and see
how we can solve it using gauss seidel method we
will just start with three four iteration and then I
will tell the numerical scheme for a one d unsteady conduction problem and the rudiments
of finite volume method also I will explain in
the next class.
