welcome back to the last league of module 3 on finite difference method we are working

on the applications of finite difference method to scalar transport and its computer applications

now let us have a recap of previous lecture the lecture 8 which we had earlier we completed

self-centered grid development for finite difference discretization for 1-d heat conduction

and we also discussed how this simple approach can be easily extended for multi dimensional

problems

then we discussed application of finite difference method to 2 dimensional heat conduction transient

heat conduction and today we are going to focus in todays lecture on advection diffusion

problem and computer implementation aspects

the brief outline of what we will do today application of finite difference method to

1 dimensional advection diffusion problem extension to multi dimensional should be fairly

simple we will have a brief look at different schemes which are available different difference

scheme for advection terms and we will then focus on the computer implementation of the

algorithms which we learnt so far for real life cfd

primarily we will have a look at the choice of languages which are used programming paradigm

and data structures

now let us first take 1 dimensional advection diffusion case so if we are dealing with steady

state problem 1 dimensional advection diffusion for a scalar variable phi is governed by d/dx

of rho u phi this left hand side this term we would call as convection term because this

represents the advection or convection of a scalar phi due to velocity u on the right

side we have got d/dx within bracket gamma d phi/dx where we will presume that gamma

may be a function of the spatial variables so this term we refer to as diffusive term

now this usual symbols u stands for specified velocity rho is density and what follows we

are going to assume that these might be considered as constant and gamma is diffusive array now

let us take this typical boundary value problem suppose phi values are specified at both ends

with domain of length l then boundary conditions are phi at x=0=phi0 where phi0 is a specified

value and similarly the right end phi at x=l is phi subscript l

this phi 0 and phi l they would both denote these specified values incorporation of gradient

boundary conditions should be exactly done in the same way as we have discussed earlier

for 1 dimensional conduction problem

so let us see your choice of schemes which you can have we will primarily choose central

difference scheme for finite different approximation of diffusive term that is the second hour

derivative which we had on the right hand side for convective terms on the left hand

side we can use cds that is central difference scheme or upwind scheme and once we have done

our discretization

a discretized equation for an interior node i can be represented in our standard template

form that ap superscript i that is coefficient a for current point p or i-th node phi p+awi

phi i-1+instead of phi this should be capital a and this should be small phi phi at i+1=qi

now let us note that these coefficients which we have put ap aw and ae they contain the

contributions coming from the finite difference discretization of both convective and diffusive

terms

that is to say we can break it into 2 parts a=a superscript d+ac and we will have a look

at how do we get these components with different choices of our discretization scheme

now let us first have a look at the contribution which comes from the diffusive term

we have represented our discrete equation at node i as

an interior node node i let us draw a standard grid of the

left end x=0 and this is x=l and if we use our vertex centre by a different scheme these

nodes will be numbered 1 2 and so on so i-th node i-1 i+1 and the last node would be at

n+1 if we had n number of divisions

so we wrote our generic discrete equation at node i as ap the subscript p stands for

our grid node i to reinforce to put superscript i as well so api phi i+awi w stands for the

i-1th grid+aei phi i+1=q1 and we said any of these components at a they will consist

of 2 parts one which is contributed by the diffusive term and another one which is contributed

by the finite difference discretization of convective term

so now let us first take the case of diffusive term and how do we discretize it we have already

seen discretization in detail so we will just put a summary here

we would use central differencing so what we are looking at we are looking at a term

del/del x of gamma *del phi/del x at node i and we said look this we would discretize

using the central difference formula so we are focusing on our grid node i

so let us chose these 2 intermediate node points i+1/2 and i-1/2 so in terms of these

values the outer derivative can be approximated as gamma*del phi/del x at i+1/2-gamma*del

phi/del x at i-1/2/x of i+1/2-xi-1/2 now there are derivatives which appear in the numerator

we would again use cds cds for first order derivatives

in a numerator so that will give us gamma del phi/del x i+1/2

this can be approximated in terms of the values at node point i+1 and i because i+1/2 is the

midpoint i and i+1/2 gamma i+1/2 phi i+1-phi i let me reward back to delta notation to

delta xi

similarly del phi/del x i-1/2 this can be approximated as gamma i-1/2 phi fi-phi fi-1/

delta xi-1 we can substitute the expression for these expanded variables or the expansion

for the derivatives appearing in numerator and thereby we can obtain the final form for

the diffusive term del phi/del x gamma*del phi/del x of i gamma i+1/2 phi fi+1-phi i

delta xi-1-gamma i-1/2 phi i-phi fi-1 delta xi this whole

thing divided by 1/2

this xi+1/2 and xi-1/2 their difference can be represented in terms of delta xi and delta

xi+1/2 or rather we can simply write it as 1/2xi+1--1 delta xi-1 you can simplify this

expression further and that i would leave it as an exercise to simplify

this expression that is all that you need to do is collect terms in phi with same subscript

together and obtain expressions for apd awd and ae of d

so this i would leave as an exercise for you it is a fairly simple exercise so this apd

would basically you will obtain by collecting the terms which multiply phi i+1 and awi would

be the terms or the coefficient which multiplies after gathering the appropriate terms which

multiply phi i-1 similarly aed that would represent the term which form the coefficient

of phi i +y

next let us move on to convective term and the convective term we have many choices available

we can use cds scheme for an advective term when i just like to make a simple statement

or generic statement that it can lead to what we call spurious oscillations that even though

oscillation might be smooth function we would find lots of oscillations if the so called

cell peclet number>2

and to ameliorate the situation what we normally do is we use upwind scheme for advective term

and there are various choices available we will see one of them today that is the first

order of upwind scheme this term it eliminates the oscillations altogether but it introduces

what we call a diffusive term in the equation which gives us a highly damped profile so

that is why people say that this particular scheme is highly diffusive

what do you mean upwind scheme let us have a brief look at first cds and then the upwind

scheme

so discretization of advective term

there are various options option 1 what we did for the diffusive term we had used second

order accurate cds scheme so shall we use the same scheme for our convective term as

well so if we use the central difference approximation then what we would get the term which we are

looking was d/dx of rho u phi at node i since we are looking at node i in central difference

approximation and terms involved would be at node i+1 and i-1

so this can be approximated as rho u phi at i+1-rho u phi at i-1 divided by the spacing

between these 2 points that is x of i+1-xi-1 so what is the contribution which we get for

the convective term at apc and this case would be our contribution to the i-th node term

apc=0 because on the right side of this equation there is no term which contains phi i and

how about phi i-1-th term coefficient that is what we said that will give us the contribution

to the aw term what will that be it will definitely be non 0

so can you rearrange and find the way out this i would again leave as a very simple

exercise so find out what will be the contributions to this apc and awc and aec if you use central

difference approximation now the remark which we had earlier made was node though the cds

is second order accurate that is accuracy is second order but it can produce oscillatory

results

if local we also use sometimes a term cell local peclet number how do we define this

local peclet number use a symbol pe this is defined as rho u delta x/gamma so this is

definition of local peclet number if this is greater than 2 we are in trouble when we

use our cds scheme so what is the solution then to get rid of oscillations rid of spurious

numerical oscillations

i want to reinforce the term spurious numerical because such oscillations are not actually

the part of our exact solution of the problem they come here simply because of our choice

of the numerical approximation or finite difference approximation we have to use what we call

use upwind difference schemes

now here we have used the term upwind what do you mean by upwind upwind comes from normal

connotation of when we say the wind is blowing so this is an advection diffusion problem

and advection depends on the velocity so i suppose this focuses on a specific point in

a space it could be a grid point i look at the velocity here now there are 2 possibilities

so x direction that u could be >0 or 0 would mean what that is our velocity but

flow is going from left to right

so what would be our upwind direction upwind means that we go against the wind or

against the flow velocity so upwind direction would be

towards negative x direction so this is one case so now in this case upwind difference

scheme should involve what this would involve values at the nodes which are not towards

the right but in a state they are towards the upwind or upstream side

so some people also use in place of upwind they also use this term very often upstream

both are synonymously used so in this case if u is >0 that is flow is in positive x direction

upwind differencing should involve values at upstream nodes

that is which are the upstream nodes they will also involve the values of node i that

is i-1 i-2 and so on along with values at node i

so these are the nodes which we will have in this case which should be involved in our

formulation and the physical logic given is that is the effect of the velocity would be

felt more because upstream whatever happens at the upstream nodes that affects the things

happening at node i so now our derivative approximation should involve the values at

these nodes which are towards the left hand side in this case

now there could be another possibility the second case would be that u is<0 that is flow

is directed in negative x direction so now in this case which ones will be our upwind

nodes this is the direction of the flow at a given node location i so upstream nodes

would be nodes i+1 i+2 i+3 and so on so now in this case our upstream nodes

or upwind nodes are i+1 i+2 i+3 and so on

so use their values function values at these nodes in finite difference approximation

now what are the choices which are available to us in the case of upwind difference schemes

the simplest choice is use of a single upstream node first order

upwind schemes which would essentially correspond to our first order what we call forward difference

or backward difference schemes so here basically we will approximate that d/dx at rho u phi

that should be approximated as what if our u is > 0 then we will have to take the values

to the left of i that is we will have to go for backward differencing

we will have rho u phi at i-rho u phi at i-1/xi-x of i-1/ and if u were >0 that is our flow

is directed from right to left so in that case we have to consider our upwind nodes

becomes i+i-th node so we will have the values rho u phi at i+1-rho u phi at i/x of i+1-xi

so in the first case we have used what we call at backward difference approximation

and when used <0 we would have used what is exactly our forward difference scheme

now what are the nice properties of the scheme no oscillations so this advantageous feature

no oscillations in numerical solutions for any value of pe it does not matter how

smaller how larger peclet number that is how strong is the convective term our first order

upwind scheme would always give us a smooth solution and downside is the disadvantage

which we should be aware

that is excessive numerical diffusion or some people also prefer

to call it a damping though it is very difficult to say that why is this particular term damping

is used in this case

and just to give you a view of graphical illustration suppose we are dealing with the advection

diffusion of the sharp front this is phi along this we have got x the sharp front is supposed

holds to get a constant value up to certain point and suddenly

takes a jump and the values go to 0 now if you use so this is our so called exact solution

if we use a cds based scheme it could give us this sort of a solution lots of oscillations

but if we use upwind difference scheme we might have a fairly smoothed out solution

so you can easily say that this is as if we have introduced additional amount of diffusivity

in our system now is there any way in which we can improve the situation one-way improvement

which you can immediately think of use a higher order upwind scheme so improvements higher

order upwind schemes for example a second order accurate fds or bds that is forward

difference or backward difference scheme

but surprisingly this is also prone to oscillations so this is not real improvement then our next

thing is to go for developed recently which are called high resolution

tvd schemes so this work similar to our upwind scheme but we are not going to discuss these

in detail at the moment i would in fact leave it as reading assignment to you that is what

is a summary that first order upward scheme it works beautifully in damping out or in

eliminating this spurious oscillation but it is highly diffusive

second order upwind scheme based on a one sided fds or bds it is again prone to oscillatory

behavior so in usual or commercial cfd course we go for what we call total variation diminishing

or tvd schemes this one possibility and other set of schemes are called essentially non-oscillatory

schemes and our first order upwind scheme in fact is a member of these generic tvd schemes

so for details i would suggest that you refer to the books by chung and versteeg and malalasekera

on computational fluid dynamics i will give you the full reference of these books towards

end of the lecture

now let us move on to our computer implementation aspects so if you want to develop a cfd code

we have aspects which we must consider that is at first which algorithms we are going

to chose to solve our partial differential equation governing a specific cfd problem

so we have primarily these 3 popular choices finite difference method finite volume method

and finite element we have already discussed in this module in detail finite difference

method

so in the next lecture we will discuss development of a simple one dimensional code based on

finite difference so in set of algorithm first choice is our discretization method and the

next is the solvers for discrete system so we have to make these choices in fact in our

useful cfd code we have got to provide for a set of solvers for the discrete systems

then which language to choose and which programming paradigm

the most popular languages in cfd applications have been fortran and fortran has been used

since beginning since 1960s and more recently c and c++ have been claiming the space which

was earlier exclusively occupied by fortran then should we go for procedural approach

which is the way the fortran works or should be go for object oriented code which gives

us what we call as easily modifiable program

so we have to make a choice of these 2 and then the last thing which we have to think

of is the choice of the data structures should we go for fixed size arrays or should we go

for the simple records arrays or we go for what we call some more complicated classes

structures or records class and struct they are part of a c language records are now available

as a feature in modern fortran language

so should we go for these ones bit more complicated ones or we go for simple arrays the next one

should we go for the fixed sized arrays that is to say we have when we are writing the

program we would think of a size or maximum number of nodes which 1 might use in discretization

make all arrays of that size or we can allow the array to grow that is for different problem

different users might chose different grid sizes and then the grid size can be input

and based on the grid size the sizes of these arrays the structures or records could be

considered

now as far as the algorithms are concerned whether you are going to chose finite difference

finite volume or finite element we will discuss finite volume and finite element in bit more

detail in later modules it would depend on your preference majority of the commercial

codes in cfd are based on either finite difference or finite volume method there are few codes

available on fem as well but it is your call

if you are developing your own code which one you are more comfortable with and this

is for the common practice today that fdm and fvm finite difference and finite volume

method they are more popular in commercial cfd applications locations and finite element

is preferred in structural analysis there is no reason why one cannot go for fem for

writing a cfd code now next is solvers for discrete algebraic systems

in our usual code we should provide a set of iterative solvers and of course these iterative

solvers require what we call pre conditioners now these 2 aspects we will consider in our

next module on solution of algebraic equations now i have not mentioned any direct solver

because in our usual cfd when the grid size are pretty large number of nodes are very

large we would rarely if ever use a direct solver

similarly if you go to time dependent problem we have to provide for a set of one step and

multi step time integration methods both of implicit types and explicit types so that

depending on the problem the user can opt for either implicit one or an explicit approach

now i would just like to point out here that as it regards the choice of their language

you can chose either c c++ or fortran but please remember if you want to become a cfd

developer a basic knowledge of both these languages is required the reason is very simple

fortran has been used in fact fortran language was developed exclusively for scientific programming

and so lots and lots of work is available in fortran

there are many libraries which are available in fortran majority of the codes which are

available in commercial and free domain they are based on fortran and if you want to make

use of those libraries the user must know at least must have a working knowledge of

fortran and a good knowledge of c or c++ now regarding procedural object oriented programming

i would like to make one small comment that object oriented code is what we call it provides

for better software engineering in the sense that the codes were easily modifiable easily

extensible because you will work in terms of classes

nuclear classes could be derived based on those classes without affecting the current

implementation so that is why it is called a better software engineering but this object

orientation has got severe performance overheads and cfd applications are what we call highly

compute intensive so we would like to minimize these non essential overheads to the extent

possible

hence in cfd programming our prefer programming paradigm is procedural programming we would

not use in fact majority of the software developers in cfd they do not use object oriented code

procedural programming is preferred now few more guidelines or what i would say the basic

guidelines for software development whichever language you chose whether you chose

fortran or c and c++

there are few basic thumb rules which you should follow so basic guidelines for development

of cfd code in fact these guidelines apply they are taken from software engineering so

they apply for development any computer code and they are equally pertinent to our cfd

applications so the first approach which we should use is what we call design

your code as we call it modular

that is to say the cfd code consists of a set of modules which can be independently

developed

and tested before combining together

this is very similar to a normal engineering practice for instance if you want to develop

or we want to design a car or manufacture a car what do we do the entire car is not

manufactured at one go different subsystems are designed separately they are manufactured

separately tested separately and once each component has been tested then it is assembled

together to get our final product that is our car

similarly if we want to develop a nice cfd code we should break our big code into small

small modules or pieces develop each module separately test it and only then combine it

together and second guideline which you should follow is document

your code now by documentation there should be 2 set of documentations

the first one is your design document which will contain a description of algorithm

and pseudo code

in hard copy form

and the second and which is equally important is provide liberal comments comments or set

of statements or set of lines in a code which are there only for human programmer they are

ignored by the complier so we have to provide liberal comments in our code that is specifically

for each routine or function in fortran we use the terms routine or sub routine in c

or c++ we call them as functions for each sub routine or functions one should provide

in our first what that function does objective of the function second description of

input and output parameters and you should also put a date stamp

which would help you to keep track of when that particular routine was return or modified

so the first guidelines that you should develop your code as a modular code and the second

aspect which we had was a documentation

the third aspect which is equally important is incremental development and testing

so there are 2 ways you are very ambitious you know the knowledge of your domain you

thoroughly know your computer programming language you also know your cfd algorithms

and you might be tempted to go for write a big piece of code and then try and run it

it might in all probability 99% chances would be that your code would fail and you would

find it very difficult to detect where you had made an error

so what you should do is we should take an incremental step by step approach and its

incremental which we can also call a step by step that is we add one routine at a time

add one routine or function at a time test this function for a most possible or most

likely values of parameters to verify that whatever algorithm has been coded in that

function that works correctly and secondly for extreme values

the parameters to detect which are the likely situation and whichever algorithm might fail

and what sort of error statements or exact cases we need to incorporate in our code so

if you follow this incremental approach you can starting from a code of let say 10 lines

you can easily built it step by step into a code which might run into millions or billions

of life so even your very big commercial code which might contain millions of lines that

has been built up using this incremental development and testing approach

there might have been many program who are involved in development but each programmer

was assigned task of developing a set of routines or set of functions or set of modules and

in developing their modules or function they had followed or they had instructed in this

industry to follow this incremental developmental and testing approach the code has to be written

in a structured format it must be documented properly both in hard copy form as well as

inside the code

and then unit testing it ensures that our code has been debugged properly and it is

ready for what we call final integration so once we have unit tested thoroughly the task

of integration would become a lot easier and we will get a reasonably error free code we

cannot say that it is a code which have developed or product which we have got that is perfect

that is 100% error free that can never happen in a human endeavor but by following this

incremental approach we can minimize the occurrence of those error

so we can confidently say that we have developed a reasonably error free code we are going

to follow these guidelines and in developing a very simple one dimensional heat conduction

code based on finite difference method that is what we are going to discuss in the next

lecture we will primarily focus on our approach so that you can use the same approach to develop

bigger codes let say multi dimensional problems

the references which we have refer today specifically for some algorithms there are 2 which you

can refer to computational fluid dynamic second edition cambridge university press by chung

tj which is published in 2010 similarly you can also refer to ferziger and peric this

book is published in 2003 computational method fluid dynamics in the next lecture towards

end i will give you consolidated set of the books or references which would be very useful

for learning finite difference method and the cfd application based in that