P S Ghoshdastidar

Associate Professor
Department of Mechanical Engineering
lIT, Kanpur
The present book covers the fundamentals of what is commonly known as
Computational Fluid Oynamics (CFO). The past two decades have witnessed a
phenomenal growth in this area due to the developments in the field of computers.
CFO has now become an integral part of the engineering design and analysis.
Engineers can make use of the CFO tools to simulate fluid flow and heat transfer
phenomena in a design and predict the system performance before manufacturing.
The advantages of CFO are numerous, namely, fewer iterations to the final
design, shorter time to launch the product, fewer expensive prototypes and so on.
Furthermore, CPO provides a cost-efficient means of testing new designs and
concepts that would otherwise be too expensive and hazardous to investigate.
Much of the material in this textbook has been used in a post -graduate course at
the Indian Institute of Technology , Kanpur for over a decade. It is assumed that the
reader has an adequate undergraduate background in Heat Transfer, Fluid Flow,
Calculus and Computer Programming in FORTRAN. The book is suitable as a
text for a one-semester course at the post-graduate or advanced undergraduate
level. It can also be used for self-study by practising engineers.
The book primarily follows finite difference method of discretization.
However, in the Appendix A, other important schemes such as finite element and
finite volume are also discussed. An emphasis has been laid on the physical
understanding of the problems. Most of the methods have been illustrated with
detailed example problems and the solution procedure. Several exercise problems
are given at the end of various chapters. Readers are encouraged to solve these
problems, to get a better understanding of various numerical techniques discussed
in the book. Chapter 7 gives details of two new numerical methods and their
applications. Chapter 8 illustrates the application of CFO in solving industrial
problems. An important subroutine (TOMA) in which tridiagonal matrix
algorithm is programmed is listed in Appendix C.
The softcover version of this book also contains a floppy diskette. The diskette
contains 21 files comprising 10 programs, I subroutine and 10 output files. The
programs are given in FORTRAN language and can be run on a PC-AT or
Pentium as well as on mainframe computer systems having a FORTRAN 77
compiler. Basically, the floppy contains programs and solutions to some unsolved
viii Preface

problems in the book and solutions '[programs] to some typical problems

discussed in this book.
I would like to acknowledge the interaction with the students, both in and
outside the class, which has greatly contributed towards the shaping of this book.
Their suggestions and comments have been useful in writing this text. I have also
been benefitted by the lively discussions with some of my colleagues. Special
thanks are due to the post-graduate students Suresh Singh, Vipin Kumar,
R Mahesh Kumar and Kali Sanjay who have assisted me in developing some of
the computer programs in the floppy diskette.
CONTENTS
I also wish to acknowledge the support and encouragement provided by the
editorial and production team of Tata McGraw-Hill, particularly, Ms Vibha
Mahajan and other members of their highly skilled editorial team. I am also grate-
ful to the anonymous reviewer whose valuable comments and suggestions for
Preface vii
improvement have gone a long way in the formation of the final version of this
book.
The typing was carried out with great care and patience by U S Mishra. The
1.mTRoDucnoN 1
figures were drawn with great competence by B N Srivastava. 'J.ne writing of this 1.1 What is Computer Simulation? 1
book would not have been possible without the generous financial support ofthe 1.2 Advantages of Computer Simulation 2
Curriculum Development Cell under the Quality Improvement Programme at the 1.3 Applications of Fluid Flow and Heat Transfer 2
Indian Institute of Technology, Kanpur. 1.4 Why is Computer Simulation Necessary in Fluid Flow
Last but not the least, the greatest contribution to this work has been the and Heat Transfer? 3
patience and encouragement of my wife Sumita and my daughter Shreya who 1.5 Basic Approach in Solving a Problem by Numerical Method 4
often withstood my moody behaviour during the writing cf this book with smiling 1.6 Problem Complexity 4
faces. Gratitude not expressible in words is due to my parents for their blessings 1.7 A Comparative Study of Experimental, Analytical
and good wishes. and Numerical Methods 5
1.8 Methods of Discretization 7
1.9 Justifications for the Choice of the Finite Difference
P S GHOSHDASTIDAR Method 8
Summary 8
References 9

2. PARTIAL DIFFERENTIAL EQUATIONS 10

2.1 Classification of PDEs 10
2.2 Elliptic, Parabolic and Hyperbolic Equations 11
2.3 Initial and Boundary Conditions 13
2.4 Initial and Boundary Value Problems 14
2.5 Conventions Followed in This Text 15
2.6 What About Hyperbolic Problems? 15
2.7 How Many Initial and Boundary Conditions do You
Need for Completely Defining a Problem? 16
Summary 17
References 17
Exercise Problems 17
x Contents Contents xi

3. INTRODUCTION TO FINITE DIFFERENCE, 4.22 False Transient Approach 89

NUMERICAL ERRORS AND ACCURACY 21 4.23 Three-Dimensional Transient Problems 89
4.24 Problems in Cylindrical Geometry: Handling
3.1 Introduction 21 of Condition at the Centre 90
3.2 Central, Forward and Backward Difference 4.25 Problems in Spherical Geometry 95
Expressions for a Uniform Grid 21 4.26 One-Dimensional Transient Conduction in
3.3 Central Difference Expressions for a Non-Uniform Grid 25 Composite Media 96
3.4 Numerical Errors 28 4.27 Treatment of Nonlinearities in Heat Conduction 98
3.5 Accuracy of Solution: Optimum Step Size 29 4.28 Irregular Geometry 106
3.6 Method of Choosing Optimum Step Size:
Summary 110
Grid Independence Test 29
References 111
Summary 30 Exercise Problems 112
References 30
Exercise Problems 30 5. NUMERICAL METHODS FOR INCOMPRESSmLE
FLUID FLOW 119
4. NUMERICAL METHODS FOR CONDUCTION
HEAT TRANSFER 32 5.1 Introduction 119
5.2 Governing Equations 120
4.1 Applications of Heat Conduction 32 5.3 Difficulties in Solving Navier-Stokes Equations 120
4.2 Steady and Unsteady Conduction 32 5.4 Stream Function-Vorticity Method 121
4.3 Dimensionality in Conduction 33 5.5 General Algorithm for Solution by If! - ~ method 123
4.4 Some Important Examples of Heat Generation in a Body 33 5.6 Creeping Flow (Very Small Reynolds Number) 125
4.5 How Does the Classification of a Conduction 5.7 InviscidFlow(Steady) 127
Problem Help? 33 5.8 Determination of Pressure for Viscous Flow 128
4.6 Basic Approach in Numerical Heat Conduction 33 5.9 Is the Transient Approach Used for Solving Steady
4.7 One-Dimensional Steady State Problem 33 Flow Problems? 131
4.8 Two-Dimensional Steady State Problem 54 5.10 If! - ~ Method for 3-D Problems 131
4.9 Three-Dimensional Problems 62 5.11 The Primitive- Variables Approach 131
4.10 Transient One-Dimensional Problem 63 5.12 Simple (Semi-Implicit Method for Pressure-Linked
4.11 Accuracy of Euler, Crank-Nicholson and Pure Equations) Procedure ofPatankar and Spalding (1972) 132
Implicit Method 70 5.13 Computation of Boundary Layer Flow 138
4.12 Stability: Numerically Induced Oscillations 71 5.14 Similarity Solutions of Boundary Layer Equations:
4.13 Convection Boundary Condition 76 Shooting Technique 140
4.14 Stability Limit of the Euler Method from Physical 5.15 Finite Difference Approach 144
Standpoint 77 5.16 Von Mises Transformation 147
4.15 Mathematical Representation of All Three Methods by
Summary 148
a Single Discretization Equation 78
References 149
4.16 Physical Representation of All Three Methods 79
Exercise Problems 149
4.17 Advantages and Disadvantages of Each of the
Three Methods 79
) 6. NUMERICAL METHODS FOR CONVECTION
4.18 How to Choose a Particular Method 80
HEAT TRANSFER 153
4.19 Consistency 80
4.20 Two-Dimensional Transient Problems 81 6.1 Introduction 153
4.21 Example Problem: Two-Dimensional Transient 6.2 Convection-Diffusion (Steady, One-Dimensional) 154
Heat Conduction in a Square Plate 83 6.3 Convection-Diffusion (Unsteady, One-Dimensional) 158
xii Contents
Contents xiii
6.4 Convection-Diffusion (Unsteady, Two-Dimensional) 159
A.7 Finite Element Method Based on Galerkin's Weighted
6.~ Computation of Thermal Boundary Layer Flows (Part A) 162
Residual Approach 268
6.6 Computation of Thermal Boundary Layer Flows (Part B) 163
A.8 Extension to 2-D and 3-D Problems 278
6.7 Transient Free Convection from a Heated Vertical Plate 170
A.9 Accuracy of Solution 279
6.8 Problem Statement 171
A.I0 Summary 279
6.9 Free Convection in Enclosure 174
A.ll Control Volume Formulation (also Known as Finite
Summary 175 Volume Method) 280
References 175 A.12 Closure 287
Exercise Problems 176
References 287
7. SPECIALTOPICS I NEWMEmODS 179 Appendix B Runge-Kutta Methods 288
7.1 Introduction 179 B.l The Essence of Runge-Kutta Methods 288
7.2 New Method (I): Application of ADI Method to Solve B.2 Simultaneous Ordinary Differential Equations 289
the Problem of 2-D Transient Heat Transfer from a B.3 Solution of Higher Order ODE by R-K Methods 290
Straight Composite Fin 179 References 290
7.3 New Method (II): Alternative to Upwind Scheme-
Application of Operator-Splitting Algorithm to Solve Appendix C Listing of Subroutine IDMA 291
Convection- Diffusion Problems 183 C.l Subroutine TDMA 291
Summary 189 C.2 Demonstration Program Showing Application
References 189 of Subroutine TDMA 292
Reference 292
8. SPECIALTOPICSn 190
Appendix D Numerical Method for Radiation in
8.1 Introduction 190
Enclosure with Diffuse-Gray Surfaces:
8.2 Transient Combined Mixed Convection and
Radiation from a Vertical Aluminium Fin
The Absorption Factor Method 293
(Ghoshdastidar and Raju, 1992) 190 . Reference 295
8.3 Heat Transfer in Rotary Kiln Reactors 205 Index 296
8.4 Modelling of Thermal Transport Processes in
Single-Screw Plasticating Extruder with Applications to
Polymer and Food Processing 220
8.5 Heat Transfer in Metal and Alloy Solidification 237
8.6 Cooling of Electronic Equipments 246
Summary 252
References 252
Appendix A Other Discretization Methods 255
A.l The Essence of the Finite Element Methods 255
A.2 Finite Element Method Based on Variational Calculus 256
A.3 Convection Boundary Condition 265.
A.4 Two-Dimensional Steady State Problem 266
A.5 Three-Dimensional Problems 267
A.6 Summary 268
 1
INTRODUCTION

1.1 WHAT IS COMPUTER SIMULATION?

The simulation of an industrial system on computer involves mathematical

representation of the physical processes undergone by the various components of
the system, by a set of equations (usually differential equations) transformed to
difference equations which are in turn solved as a set of simultaneous algebraic
equations.
At this stage, a reader uninitiated into the numerical methods may ask the
question "What is the role of computer here?". The aforesaid query is a valid one.
Many seem to forget that some of the numerical schemes (e.g. Finite-Difference)
that are extensively used today for solution of problems on computer were devel-
oped when computer was not even invented. Now, to return to the original ques-
tion, the answer is that with the aid of the algorithm of the solution method trans-
lated into a programming language like FORTRAN fed into a computer which
does the arithmetic operations at a tremendous speed one can obtain the solution
of mathematical equations in seconds or even in fraction of a second. A simple
example will clarify this point. One can very easily solve a set of three linear
simultaneous algebraic equations. by hand through Gaussian Elimination
method.Typically in this method, for a system ofn equations the total number of
multiplications and divisions is roughly.!. n3• Therefore, for n = 3, the number of
3
operations is 9, which is clearly manageable by hand calculations. However, for
= =
n 10, this number jumps to 333. For n 100, the number skyrockets to 333000.
A mainframe computer (e.g. VAX 8810) having an average megaflop * rating of
around 1 (that is, 106 arithmetic operations per sec) will solve the aforesaid prob-
lem in 0.333 second. A personal computer (e.g., IBM PC) with a megaflop rating

* A mega is a million and flops is an abbreviation for floating point operations per
second. A floating point operation is an arithmetic operation (addition, subtraction,
multiplication and division) on operands which are real numbers with fractional
parts. NormaUy multiplications and divisions are counted asinajor arithmetic op-
erations as compared to addition and subtraction on a computer.
2 Computer Simulation of Flow and Heat Transfer Introduction 3

of 0.01 will take about 33 seconds. In computer simulation, it is possible to han- • Design of IC engines.
dle one hundred or more (even greater than thousand) number of equations and • Optimization of heat transfer from cooling fins.
there lies the necessity of using the computer. • Aircraft and spacecraft.
• Design of electrical machinery and electronic circuits.
1.2 ADVANTAGES OF COMPUTER SIMUlATION • Cooling of computers.
• Weather prediction and environmental pollution.
Now that the use of computer in simulation is established, let us enumerate some • Materials processing such as solidification and melting, metal cutting,
of the important advantages of computer simulation (also known as numerical welding, rolling, extrusion, plastics and food processing in screw
simulation): extruders, laser cutting of materials.
• It is possible to see simultaneously the effect of various parameters and • Oil exploration.
variables on the behaviour of the system since the speed of computing is • Production of chemicals such as cement and aluminium oxide.
very high. To study the same in an experimental setup is not only difficult • Drying.
and time-consuming but in many cases, may be impossible. • Processing of solid and liquid wastes.
• It is much cheaper than setting up big experiments or building prototypes • Bio-heat transfer (as in human and animal bodies).
of physical systems. It is no wonder that J.B. Joseph Fourier, father of the theory of heat diffusion
• Numerical modelling is versatile. A large variety of problems with made this remark in 1824 "Heat, like gravity, penetrates every substance of the
different levels of complexity can be simulated on a computer. universe; its rays occupy all parts of space. The theory of heat will hereafter form
• "Numerical experimentation" (another synonym for computer simulation) one of the most important branches of general physics". Lord Kelvin, in 1864
allows models and hence physical understanding of the problem to be obtained a rough estimate of the age of earth based on an idea proposed by Fourier
improved. It is similar to conducting experiments. in 1820 to be 94 x 106 years (0.094 billion years) by applying the principle of
• In some cases, it is the only feasible substitute for ~xperiments, for transient heat conduction in a semi-infinite solid. Modem dating methods have
example, modelling loss of coolant accident (LOCA) in nuclear reactors, revealed the age of earth to be approximately 4.7 billion years. So Kelvin's result
numerical simulation of spread of fire in a building and modelling of was not really too far off the mark considering the fact that the data for the meas-
incineration of hazardous waste. ured value of the geothermal gradient (rate of increase of temperature of earth
However, it is to be emphasized that not every problem can be solved by with depth), average thermal diffusivity of rock and the initial temperature of
computer simulation. Experiments are still required to get an insight into the molten earth when cooling began available with him at that time were not very
phenomena that are not well understood (and hence cannot be translated into the ,ccurate. The aforesaid example is probably the first known application of heat
language of mathematics) and also to check the validity of the results of computer transfer simulation.
simulation of complex problems.
1.4 WHY IS COMPUTER SIMUlATION NECESSARY
1.3 APPUCATIONS OF FLUID FLOW AND HEAT IN FLUID FLOW AND HEAT TRANSFER?
TRANSFER
Ifope looks at a classical textbook on fluid dynamics and heat transfer, one would
Fluid flow and heat transfer playa very important role in nature, living organisms find only a handful of analytical (or exact) solutions. In actual situations, prob-
and in a variety of practical situations. More often than not, flow and heat transfer lems are lot more complex as in those involving non-linear governing equations
are coupled and rarely an engineer solves a problem of either pure fluid flow or and/or boundary conditions, and irregular geometry which do not allow analyti-
pure heat transfer. In many applications flow and heat transfer are accompanied cal solutions to be obtained. Therefore, it is necessary to use numerical tech-
by chemical reaction and/or mass transfer. The various applications of fluid flow niques for most problems of practical interest. Furthermore, to design and
and heat transfer are: optimize thermal processes and systems, numerical simulation of the relevant
transport phenomena is a must, since experimentation is usually too involved and
• All methods of power production, e.g. thermal, nuclear, hydraulic, wind,
expensive. However, necessary experimentation must still be done in checking
and solar power plants. the accuracy and validity of numerical results. Sometimes, numerical model can
• Heating and air-conditioning of buildings. be refined by input from results of a companion experimental set-up for the same
• Chemical and metallurgical industries, e.g. furnaces, heat exchangers, problem.
condensers and reactors.
4 Computer Simulation of Flow and Heat Transfer Introduction 5
1.5 BASIC APPROACH IN SOLVING A PROBLEM For the problem shown in Fig. 1.1, since there are 16 tiny squares, GG = 16. If
BY NUMERICAL METIIOD a direct method like Gaussian Elimination method is used, then SS = 1, VG = 1,

Suppose we wish to. obtain the temperature field in the domain as shown in FP = 1. (9)3 = 243.
3
Fig. 1.1. We imagine that the domain is filled by a grid, and seek the values of
temperatures at the grid points. Therefore, the energy equation (which is the PC = 16 x 1 x 1 x 243 = 3.89 x 103 operations
governing differential equation for the problem) is valid at all the grid points. The Thus, even with a fairly coarse grid, the number of operations is quite large.
governing differential equation is then transformed into a system of difference Currently a large number of realistic applications like modelling of supersonic
equations resulting in a set of simultaneous algebraic equations which means that aircraft or weather prediction requires 1012 to 1014 operations per solution
if there are 100 grid points (where variables are not known) there will be 100 (Rajaraman, 1993). If each solution has to be done in about an hour, then the
equations to solve per variable. The simplification inherent in the use of algebraic average speed of computing should be 1014/60 x 60 operations per second which
equations rather than differential equations is what makes numerical methods so is equal to 27,700 Megaflops! The peak megaflop rating of modem supercom-
powerful and widely applicable. puters is around 1000 megaflops (Rajaraman, 1993). From the aforesaid exam-
x Grid points ple, it is clear why we need supercomputers with speeds in thousands of megaflops
(known temperatures) range to solve extremely complex problems.

o Grid points 1.7 A COMPARATIVE STUDY OF EXPERIMENTAL,

(unknown temperatures)
ANALYTICAL AND NUMERICAL METIIODS

(a) Experimental Method Experimental methods are used to obtain reliable

informations about physical processes that are not well understood, e.g. combus-
tion and turbulence. It may involve full scale, small scale or blown-up scale
models. The major disadvantages of experimental investigations are high cost,
measurement difficulties and probe errors. Often, small scale models do not al-
ways simulate all the features of the full scale set-up. The advantage is that it is
most realistic.
(b) Analytical Method Analytical methods or methods of classical mathemat-
ics are used to obtain the solution of a mathematical model consisting of a set of
differential equations which represent a physical process within the limit of as-
sumptions made. Only a handful of analytical solutions are available in heat trans-
Fig. 1.1 4 x 4 grid for computation of two-dimensional fer and fluid mechanics because analytical methods are inadequate in handling
steady-state heat conduction in a square plate complex boundary and non-linearities in the differential equations and/or bound-
ary conditions. Furthermore, the analytical solutions often contain infinite series,
1.6 PROBLEM COMPLEXI1Y -special functions, transcendental equations for eigenvalues, etc. the numerical
evaluation of which becomes quite cumbersome.
In general, the problem complexity is described by the formula
(c) Numerical Method As explained in Sec. 1.5, a numerical prediction works
PC= GGx VGxSSxFP out the consequences of a mathematical_ model but the solution is obtained for
where, =
PC problem complexity variables at discrete grid points in the computational domain in contrast with
GG = geometry of the grid system analytical method which gives closed form solution at all points (theoretically
VG = variables per grid point infinite number of points in the solution domain or continuum). The major advan-
SS = number of steps per simulation for solving problem tages of numerical solution are its abilities to handle complex geometry and non-
FP = number of floating point operations per variable. linearities in the governing equation and/or boundary conditions. Other advan-
tages of numerical method are briefly described below:
6 Computer Simulation of Flow and Heat Transfer Introduction 7
• Low Cost While prices of most items are increasing, cost of computing
1.8 METHODS OF DISCRETIZATION
(mainframe, mini and PC's) is going down every year.
• High Speed In the past 20 years there has been a thousand fold increase
There are several methods of discretizing a given differential equation. They are
in the speed of arithmetic operations of computers.
briefly described below.
• Complete Information A Computer simulation gives detailed and
complete information of all the variables over the computational domain. (a) Finite-Difference Method .The ilsual procedure for deriving finite-differ-
• Ability to simulate realistic conditions For a computer program, there ence equations consists of approximating derivatives in the differential equation
is absolutely no problem in having an area of size 10 km x 10 km or via a truncated Taylor series. The method includes the assumption that the varia-
10-6 m x 10-6 m, 5000 °C or -50°C, hazardous or flammable material. tion of the unknown to be computed is somewhat like a polynomial in x, y or z so
• Ability to simulate ideal conditions There is no problem in idealizations that higher derivatives are unimportant. The great popularity of finite-difference
like two-dimensionality, insulated or isothermal boundary, infinite
methods is mainly due to their straight-forwardness and relative simplicity by
reaction rate. On the other hand, it is extremely difficult, if not impossible
which a newcomer in the field is able to obtain solutions of simple problems. As
to set up the same in experiments.
a matter of fact, the subject of Computational Fluid Dynamics (commonly abbre-
The disdvantages of numerical predictions are the associated truncation error viated as CFD) was born as early as 1933 (remember that world's first computer
and round-off error and the difficulty in simulating complicated boundary condi-
called ENIAC was built at the University of Pennsylvania, USA, in 1946) with
tions.
the remarkable work (published in the proceedings of the Royal Society) ofThom
The aforesaid discussion can now be represented in a capsuled form in
Table 1.1. (1933) who solved the Navier-Stokes equations for the steady, incompressible
Table 1.1 Comparison of Experimental, Analytical and viscous flow around a circular cylinder by finite-difference method using hand
Numerical Methods of Solution calculations. However, several shortcomings and limitations of finite-difference
method came to light when researchers tried to solve problems with increasing
Name of the Method Advantages Disadvantages degree of physical complexity such as, for example, flows at higher Reynolds,
1. Experimental • Capable of • Equipment numbers, flows around arbitrarily shaped bodies, strongly time-dependent flows,
Being most required etc. (Fasel, 1978). This led to a search for and development of superior methods,
realistic • Scaling particularly in the areas where difference methods seemed to have disadvantages.
problem
These methods can be divided into two main categories (i) finite-element methods
• Measurement
difficulties and (ii) spectral methods.
• Probe errors
• High operating costs (b) Finite-Element Method (FEM) Finite-element methods basically seek solu-
2. Analytical • Clean, general • Restricted tions at discrete spatial regions (called elements) by assuming that the governing
information
which is
to simple
geometry
- differential equations apply to the continuum within each element. It is based on
integral minimization principle and provides piecewise (or regional) approxima-
usually in and physics
tions to the governing equations.
formulaform • Usually
restricted Finite-element methods were already found to be successful in solid mechan-
to linear problems ics applications. Their introduction and ready acceptance in fluid mechanics was
• Cumbersome due to relative ease by which flow problems with complicated boundary shapes
results-difficult
can be modelled, especially when compared with finite-difference methods. How-
to compute
3. Numerical • No restriction • Truncation ever, disadvantages of FEM arises from the fact that more complicated matrix
to linearity and round-off errors operations are required to solve the resulting system of equations. Furthermore,
• Ability to handle • Boundary meaningful variational formulations are difficult to obtain for high Reynolds
irregulargeometry condition number flows. Hence, variational principle-based FEM is limited to solutions of
and complicated problems
physics creeping flow and heat conduction problems.
• Low cost and Galerkin's weighted residual method, which is also another finite-element
high speed of method is a powerful method and circumvents the difficulties faced by variational
computation
8 Computer Simulation of Flow and Heat Transfer Introduction 9
calculus based FEM. Much current research is in progress in the use of this of a problem by numerical method is then given with the aid of a figure along with
method. the concept of problem complexity. A comparison of advantages and disadvan-
tages of numerical methods vis-a-vis analytical and experimental methods is then
(c) Spectral Method Spectral methods are generally much more accurate than
provided in detail. Finally, the chapter concludes with a brief description of the
simple first or second order finite-difference schemes. In spectral methods, in
various methods of discretization and the justifications for the choice of the finite
contrast with the discretization as in finite-difference methods, the approximation
difference method as the discretization scheme used in this book
is based on expansions of independent variables into finite (truncated) series of
smooth (mostly orthogonal) functions. In addition to its advantage as regards
higher accuracy, it can be easily combined with standard finite-difference meth-
REFERENCES
ods. A disadvantage of spectral methods is their relative complexity in compari-
son with standard finite-difference methods. Also the implementation of complex 1. Carslaw, H S and J C Jaeger, Conduction of Heat in Solids, 2nd Edition,
boundary conditions appears to be a frequent source of considerable difficulty. Clarendon Press, Oxford, 1959.
(d) Control Volume Formulation In this method, the calculation domain is 2. Fasel, H "Recent Developments in the Numerical Solutions of the Navier-Stokes
Equations and Hydrodynamic Stability Problems", Computational Fluid
divided into a number of non-overlapping control volumes such that there is one
Dynamics, Edited by Wolfgang Kollmann, Hemisphere, Washington, D.C., 1978.
control volume surrounding each grid point. The differential equation is integrated 3. Field, George B, Gerrit L. Verschuur, and Cyril Ponnamperuma, Cosmic
over each control volume. Piecewise profiles expressing the variation of the un- Evolution: An Introduction to Astronomy, Houghton Mifflin Company, Boston,
known between the grid points are used to evaluate the required integrals. The 1978.
result is the discretization equation containing the values of the unknown for a 4. Jaluria, Yogesh and Kenneth E. Torrance, Computational Heat Transfer,
group of grid points (Patankar, 1980). Hemisphere, Washington, D.C., 1986.
The major advantage of this method is its physical soundness. The disadvan- 5. James, M L, G M, Smith, and J C Wolford, Applied Numerical Methods for
tage is that it is not as straightforward as finite-difference method. Digital Computation with FORTRAN, Scraton, Intemational Text Book Co., 1967.
6. Kakac, Sadik and Yarnan Yener, Heat Conduction, 3rd Ed., Taylor and Francis,
1.9 JUSTIFICATIONS FOR TIlE CHOICE OF TIlE Washington, D.C., 1993.
FINITE DIFFERENCE METIIOD 7. Orszag, S A and M. Israeli, "Numerical Simulation of Viscous Incompressible
Flows", Ann. Rev. Fluid Mech., Vol. 6, 1974, p. 281.
In this book, finite-difference method is chosen as the method of discretization 8. Patankar, S V, Numerical Heat Transfer and Fluid Flow, Hemisphere,
Washington, D.C., 1980.
because (i) for a newcomer in the field of numerical simulation of flow and heat
9. Rajaraman, V, Supercomputers, Wiley Eastern, New Delhi, 1993.
transfer, this is the best method to begin with simply because of its inherent
10. Ralston, Anthony and Philip Rabinowitz, A First Course in Numerical Analysis.
straightforwardness and simplicity, and (ii) in recent years tremendous refine- 2nd Edition McGraw-Hill, New York, 1978.
ments and advances have been made by numerous researchers in the finite dif~er- 11. Thorn, A "The Flow Past Circular Cylinders at Low Speeds", Proc. Roy. Soc., A
ence method, particularly in the areas where it was known to be weak, such as 141, 1933, p. 651.
irregular boundaries or superior accuracy.
However, interested readers may look at Appendix' A' containing a brief dis-
cussion on finite element methods and control volume method. Spectral method
has been left out in the book for the sake of brevity. The readers are however,
encouraged to go through a survey paper by Orszag and Israeli (1974) which
gives a detailed treatment of the spectral method.

SUMMARY

This chapter begins with a formal definition of the term Computer Simulation
and its advantages followed by a discussion on the applications and the necessity
of computer simulation of flow and heat ~ansfer. Basic approach in the solution
16 Computer Simulation of Flow and Heat Transfer
now-a-days solution of hyperbolic heat conduction problems are receiving atten-
tion from researchers because of its application in advanced materials processing
techniques, such as, laser cutting of materials and laser surgery. Due to short time
of the pulse in all these processes, non-Fourier conduction effects are introduced
by the finite speed of propagation of heat wave in laser irradiated materials and
hence hyperbolic heat conduction equation needs to be solved to accurately repre-
sent the physical process. Other applications of hyperbolic heat conduction are
situations involving very low temperatures near absolute zero and very high tem-
perature gradient.
Hyperbolic equations of the first order arise in many gas dynamics problems
where viscosity is unimportant and the motion is unsteady. If the motion is steady
the classification of the systems depends upon the magnitude of the fluid speed,
being hyperbolic if it is supersonic and elliptic if subsonic.

2.7 HOW MANY INITIAL AND BOUNDARY

CONDmONS DO YOU NEED FOR
COMPLETELY DEFINING A PROBLEM?

The number of initial and boundary conditions in the direction of each independ-
ent variable of a problem is equal to the order of the highest derivative of the
governing differential equation in the same direction.
Consider, for example, the conduction equation written in Cartesian
coordinates for a homogeneous isotropic solid moving with velocity
 Now, the question is: When does one use a higher order difference scheme?
There is no set answer to this. It depends on the accuracy requirement of a prob-
lem and the analyst will have to use his own judgement.
The readers are encouraged to derive Eq. (3.10) and Eq. (3.11) themselves.
3.3 CENTRAL DIFFERENCE EXPRESSIONS
3.2.4 When Do You Use Forward, Backward FOR A NON-UNIFORM GRID
and Central Difference Expressions?
Non-uniform grid is used in the region where gradients of the unknown (for
• Forward difference expressions are used when data to the left of a point at example, near the wall in boundary layer flow and/or heat transfer or heat
which a derivative is desired are not available. conduction in a semi-infinite solid) are expected to be high. Therefore, by making
• Backward difference expressions are used when data to the right of the a fine grid in the region of large gradients and coarse grid in the region whether
desired point are not available. gradients are small, one can save on computational memory and execution time.
• Central difference expressions are used when data on both sides of the Figure 3.2 shows a typical non-uniform grid, the grid spacings being non-uniform
desired point are available and are more accurate than either forward or in the x-direction only.
backward difference expressions.
 Introduction to Finite Difference, Numerical Errors and Accuracy 29

(b) Truncation Error The truncation error is due to the replacement of an exact
mathematical expression by a numerical approximation. This error has already
been discussed earlier in this chapter with respect to finite difference
approximations. Basically, it is the difference between an exact expression and
the corresponding truncated form (for example, truncated Taylor series) used in
the numerical solution.
(c) Discretization Error The discretization error is the error in the overall solu-
tion that results from the truncation error assuming the round-off error to be neg-
ligible. Therefore,
Discretization error = exact solution - numerical solution with
no round-off error

3.5 ACCURACY OF SOLImON: OPTIMUM

STEP SIZE
,"
The accuracy of a numerical solution is determined by its total error which is the
sum of the round-off error and truncation error. Hence, total error = round-off
error + truncation error.
However, it is obvious that round-off error increases as total number of arith-
metic operations increases. Again, total number of arithmetic operations increases
if the step size decreases (that is, when the number of grid points increases).
Therefore, round-off error is inversely proportional to the step-size.
On the other hand, truncation error decreases as step size decreases (or as
numb,er of grid points increases).
Because of theaf-orementioned opposing effects, an optimum step size is ex-
pectedwhich will produce minimum total error in the overall solution.

3.6 MEmOD OF CHOOSING OPTIMUM STEP

3.4 NUMERICAL ERRORS - SIZE: GRID INDEPENDENCE TEST
Three most important errors that commonly occur in numerical solution are (a) A numerical analyst has to be extremely careful as regards the accuracy of his
round-off error (b) truncation error and (c) discretization error. solution. To get the most accurate solution (that is, the solution with least total
(a) Round-off Error The round-off error is introduced because of computer's error), one has to perform a grid independence test. The test is carried out by
inability to handle large number of significant digits. Typically, in single-preci- experimenting with various grid ~izes and watching how the solution changes
sion, the number of significant figures retained ranges from 7 to 16, although it with respect to the changes in grid sizes. Finally, a stage will come when chang-
may vary from one computer system to another. The round-off error arises due to ing the grid spacings will not affect the solution. In other words, the solution has
the fact that a finite number of significant digits or decimal places are retained now become independent of grid spacing. The largest value of grid spacings for
and all real numbers are rounded off by the computer. The last retained digit is which the solution is essentially independent of step sizes is chosen so that both
rounded off if the first discarded digit is equal to or greater than 5. Otherwise, it is the computational time and effort and the round-off error are minimized. A.nex-
unchanged. For an example, if five significant digits are to be kept in place,
ample of a grid independence test is given in Chapter 4.
5.37527 is rounded oftto 5.3753, and 5.37524 to 5.3752.
30 Computer Simulation of Flow and Heat Transfer

SUMMARY .

This chapter introduces readers to finite difference approximation of derivatives

of continuous and differentiable functions based on Taylor series expansion. The
central, forward and backward difference expressions for first and second deriva-
tive of a function have been derived. Difference expressions of higher order accu-
racy are also given. Next, application of non-uniform grid and a detailed deriva-
tion of schemes for first and second derivatives of function evaluated at a point in
a non-uniform grid follow. The concepts of round-off error, truncation error,
discretization error and total error have been discussed in detail with reference to
solution accuracy. Finally, method of choosing optimum grid size (grid independ-
ence test) is briefly mentioned.

REFERENCES

1. Carnahan, Brice, H A Luther, and James O. Wilkes, Applied Numerical Methods,

John Wiley & Sons, New York. 1969.
2. Constantinides, Alkis, Applied Numerical Methods with Personal Computers,
McGraw-Hill Inc., New York, 1987.
3. Jaluria, Yogesh, Computer Methods for Engineering, Allyn and Bacon, Inc.,
Boston, 1988.
4. James, M L, G M Smith, and J C Wolford, Applied Numerical Methodsfor Digital
Computation with FORTRAN, Scraton, International Text Book Co, 1967.
5. Rohsenow, Warren, M and James P Hartnett, Handbook of Heat Transfer,
McGraw-Hill Book Company, New York, 1973.
 Numerical Methods for Conduction Heat Transfer 33

4.3 DIMENSIONAll1Y IN CONDUCTION

Depending on the physical processes involved, a conduction problem may be one,

two- or three-dimensional. One and two-dimensional configurations are common.

4.4 SOME IMPORTANT EXAMPLES OF HEAT

GENERATION IN A BODY

NUMERICAL METHODS The internal heat generation may be due to flow of electrical current in a body,
chemical reaction or release of nuclear energy as in nuclear fuel rods in a nuclear
FOR CONDUCTION reactor. The source of energy may be concentrated, as in a small-diameter cur-
rent-carrying electrical conductor inside thick insulation or the source may be

HEAT TRANSFER distributed, as in the elements or shield of a nuclear reactor. The example of
uniform heat generation is the electrical heating throughout a material, as for
example in coil of an electric heater or in a microwave oven.

4.5 HOW DOES mE ClASSIFICATION OF

\~: A CONDUCTION PROBLEM HELP?
4.1 APPliCATIONS OF HEAT CONDUCTION
I The classification of a problem as mentioned earlier essentially helps the heat
Heat conduction has numerous applications in modern technology, in geological transfer analyst to decide on the approach to be used to solve the problem. It also
sciences, and in many other evolving areas such as materials processing ..Some ..tells us the level of difficulty to be encountered in obtaining the solution, that is,
examples are cooling fins or extended surfaces, solidification .and melt1O~ of .the temperature distribution within the body.
metals and alloys in metallurgical industries, welding, metal cuttl~g: quenching,
electrical chemical and nuclear heating, periodic temperature vanatlons of earth .6 BASIC APPROACH IN NUMERICAL
surface, baking oven, furnace walls, heating and cooling of buildings, to name HEAT CONDUCTION
just a few. To analyse thermal stress conditions in a mat~rial, te~perature
distribution must be obtained by solving the heat conduction equation. !he .i Many difficult problems arise in conduction, for example, variable thermal con-
starting point of all such analyses is the differential equation (energy eq.uatlOn) t
I~,

ductivity, distributed energy sources, radiation boundary conditions for which

based on the physical formulations of the phenomena relevant to conduction. ! analytical solutions are not available. Approximate solution is then obtained by
numerical method. The basic approach is to arrive at the relevant governing
4.2 STEADY AND UNSTEADY CONDUCTION differential equation based on the physics of the particular problem. They are
then converted to the required finite difference forms. To begin with, the numeri-
Since the numerical treatment of a conduction problem depends on the nature of cal solution procedure for the problem of a simple one-dimensional steady state
the conduction process, all conduction processes are divided broadly into two heat conduction in a cooling fin is described. It is to be noted that simple, closed-
categories, namely, steady and unsteady. Steady state means that temperature, form straight forward analytical solution for this problem is available. The idea is
density, etc. at all points of the conduction region is independent of time. Un- to show the use of numerical method and to compare the numerical solution with
steady state means a change with time, usually only of the te~pera~un~. Unsteady its analytical counterpart.
state problems can be further split into the categories, that IS, penodic an~ tran-
sient. Daily variation of earth's temperature due to solar eff~cts exam~hfies a 4.7 ONE-DIMENSIONAL STEADY STATE PROBLEM
typical periodic heat conduction problem. ~nother ex~mple IS the engme wall
temperature variations due to cyclic changes 10combustlon gas temperature. !he Consider the one-dimensional steady state heat conduction in an isolated rectan-
immersion of hot steel plate in a cold quenching tank is an example of transIent gular horizontal fin as shown in Fig. 4.1. The base temperature is maintained at
conduction. Transient periodic heat conduction is also not uncommon. T = To and the tip of the fin is insulated. The fin is exposed to a convective envi-
ronment (neglecting radiation heat transfer from the fin) which is at Too (Too < To).
 ·.
The above algorithm is also known as Thomas Algorithm. See Appendix C for
listing of subroutine TDMA.
Finally, it is to be noted that Eq. (4.7) might also be solved by the Gauss-
Seidel iteration scheme discussed next.
Gauss-Seidel Iterative Metbod(G-S) For a large number of equations (typically
of the order of several hundred) iterative methods, which initiate the computa-
tions with a guessed solution and iterate to the desired solution of the systems of
equations within a specified convergence criterion, using improved guesses in the
second, third iterations till the final one, are often more efficient. In this method,
unlike In direct methods like Gaussian elimination the round-off error does not
__ mulate. The round-off error after each iteration simply produces a less
IIJ LHft •••• far the next iteration. Therefore, the resulting round-off error in the
44 Computer Simulation of Flow and Heat Transfer Numerical Methods for Conduction Heat Transfer 45

This process is known as Jacobi iterative method. that is, when the system is diagonally dominant. This is also known as
Scarborough Criterion. However, convergence may still be possible even if the
Disadvantages of Jacobi Method The main disadvantage is that the computer above condition is not satisfied. Fortunately, it turns out that in fluid flow and
storage is needed for the present iteration as well as the previous one. This is heat transfer problems, finite-difference formulation indeed leads to diagonally
because all the values are computed, using previous values before any unknown dominant coefficient matrix which is the reason why for large systems Gauss-
is updated. Seidel method is so widely used.
Gauss-Seidel Method: An Improvement over Jacobi Method A significant Application of G-S Iterative Method In order to demonstrate the iteration
improvement in the rate of convergence and in the storage requirements can be process, the following system of three linear equations is solved by G-S iterative
obtained by replacing the values from the previous iteration by new ones as soon method using a pocket calculator.
as they are computed. Then, only the values of the latest iteration are stored, and
each iterative computation of the unknown employs the most recent values of the lOx, + x2 + 2x3 = 44
other unknowns. This computational scheme is known as point-by-point Gauss- 2x, + 10x2 + x3 = 51 (4.36)
Seidel method. x, + 2x2 + IOx3 = 61
 Numerical Methods for Conduction Heat Transfer 61

the second sweep, the computation again begins for the first chosen line and
sweep continues till the other boundary is reached. The same procedure is fol-
lowed for third, fourth, till final sweep when there is virtually no change in the
temperature distribution in consecutive sweeps. An alternating direction method
is also possible. In this method sweeping of the computational domian is done
alternately in x (or y) and y(or x) direction.
The sweep direction is also important in cases where in one direction one of
the boundaries has insulation or convection condition while at the other bound-
ary, Dirichlet condition is specified. For example, in the present problem
(Fig. 4.15), a right-to-left sweep would transmit the known temperature on the
right boundary into the domain; on the other hand, since no te~perature is given
on the left boundary, a left-to-right sweep would bring no such useful informa-
tion. Based on the same argument, in the y-direction, top-to-bottom sweep is de-
sirable.
The fast convergence of the line-by-line method is due to the fact that the
boundary condition information from the ends of the line is transmitted at once to
the interior of the domain in sharp contrast with point-by-point Gauss-Seidel
method in which the boundary condition information is transmitted at a rate of
one grid interval per iteration.

4.8.6 Line-by-Line Method

To alleviate the problem of slow convergence of Gauss-Seidel (G-S) method for
large number of grid points, a line-by-line method which is a convenient combi-
nation of the direct method (TDMA) for one-dimension and the G-S method can
be used. Basically, the method makes use ofthe direct nature ofTDMA and low
round-off error of G-S scheme. Fig. 4.15 Pictorial representation of the line-by-line method
In the line-by-line method, a grid line (say, in the y-direction) is chosen assum-
ing that the unknowns (say, temperature, 1) along the neighbouring lines (Le., the
4.8.7 Check for Accuracy
x-direction neighbours of the points on the chosen line) are known from their For the present problem, the accuracy of the numerical results can be checked by
latest values. Now, 1"s along the chosen line (Fig. 4.15) is solved by TDMA. comparing it with the corresponding analytical solution which is available. Sub-
This procedure is repeated for alHheJipes inrtb~:Ynd.i.re!Cti.oJl,j,J;l
the first sweep. In sequently, a grid independence test must be done to obtain the desired results.
4.12 STABILI1Y: NUMERICALLY INDUCED
OSCillATIONS

From the preceding discussion, it is apparent that all the three schemes will give
better results if time steps are made smaller. In practice, however, one would
usually like to take as large a time-step as one can to reduce the computational
effort and time. In addition to decreasing the accuracy of the solution, large time
steps can introduce some unwanted, numerically induced oscillations into the so-
lution'!11aking it physically unrealistic. Such solutions are not acceptable and the
method which produces such solution is called unstable method. This brings us to
the formal definition of a stable numerical scheme which is one for which errors
from any source (round-off, truncation, mistakes, etc.) are not permitted to grow
in the sequence of numerical procedures as the calculation proceeds from one
marching step to the next.
•
Case of One Grid Point Consider the case of only one grid point, that is, the
grid point on the insulated boundary of the plate (Fig. 4.21). Therefore,!1X. = 1,
r = ~'t', Hence, we have only one equation to solve, that is,
80 Computer Simulation of Flow and Heat Transfer

Crank-Nicholson
Advantages : More accurate than the Euler and the pure implicit methods
for the same time-step.
Stable.
Disadvantages: A set of simultaneous equations to be solved in each time-step
and hence, more computational effort.
Difficult to program on computer.

4.18 HOW TO CHOOSE A PARTICULAR METHOD

The answer is not so straightforward. Because the final choice depends on the
accuracy requirement of the problem, whether the solution at earlier or later times
is desired, and the speed and memory ofthe computer.
For example, if the solution at an earlier time is required and the accuracy can
be compromised, the Euler method should be preferred. On the other hand, if the
entire temperature-history is the one sought, then either Crank-Nicholson or pure
implicit can be used depending on the accuracy requirement.
For larger time-steps, the pure implicit scheme is better than Crank-Nicholson
as over a larger interval, the steep rise or fall in the beginning followed by a flat
tail is more realistic as compared to a linear profile throughout. For small time
interval, however, the Crank-Nicholson scheme is physically more sound
(Patankar, 1980).
Furthermore, the computer speed and memory should also be looked at. If one
is using a PC-XT, use of the Euler method won't be wise while for a mainframe
system which has much larger computer memory and higher speed, using the
Euler method will not be unwise for solving the same problem.

4.19 CONSISTENCY

A finite-difference procedure is consistent if it approximates the solution of the

PDE under study, and not the solution of some other PDE.
For example, consider the solution of a one-dime~sional transient heat conduc-
tion problem by explicit method (or Euler method).
We kne.w that the truncation error of Euler method is O[(Ll'l', (L1X)2]. Since the
truncation error tends to zero as Ll'l', M ~ 0, the explicit representation is con-
sistent with the original PDE.
On the other hand, another explicit procedure called Dufort-Frankel scheme
(which is unconditionaly stable) may not be always consistent. This can be dem-
onstrated as follows.
Dufort-Frankel scheme (Carnahan et al., 1969):
4.27 TREATMENTOF NONUNEARITIES IN
HEAT CONDUCTION
The nonlinearity in the governing differential equation or in the boundary con-
dition does not preclude its solution by one of the basic methods discussed earlier
in this chapter. One must never forget that the objective of any simple finite-
difference representation is always to approximate the nonlinear PDE by an alge-
braic equation which is linear in its unknowns. Therefore, the nonlinearities aris-
ing in the difference equation is locally linearised.

4.27.1 Nonlinear Governing Differential

Equation: Variable Thermal Conductivity
In Chapter 2, definition of a nonlinear PDE is given. In heat conduction, a
nonlinear governing differential equation results if the thermal conductivity of the
material is not a constant but a function of temperature.
Take, for example, the case of a one-dimensional transient heat conduction in
a plane wall of thickne"ssL having conductivity k = k(1) which is known. p, Care
assumed constants. The problem is pictorially described in Fig. 4.37.
104 Computer Simulation of Flow and Heat Transfer

in perfect vacuum. For example, a power plant on board a space-craft can reject
energy to the outerspace only by radiation.
However, in many problems, a medium which emits, absorbs, and transmits
radiation affects the radiative transport process. Examples of such media are
moisture and carbon dioxide. Therefore, consideration of combined convection
and radiation is important in furnaces or IC engine combustion chambers, where
the gases participate in the radiative heat transfer significantly because of the
presence of carbon dioxide, moisture and various other absorbing gases. On the
other hand, the thermal radiation from a surface to a non-absorbing and non-
emitting medium (such as air at low temperature or vacuum) or vice-versa is
independent of the convective process. Thus, the total energy transfer from or to
the heated surface, in such cases is calculated by a summation of the transfer due
to the two modes operating independently.
For low temperature levels, the radiative transfer may be neglected. While for
conduction and convection, heat transfer is almost linearly proportional to the
temperature difference, for radiation, heat transfer is proportional to the differ-
ence in the fourth power of the absolute temperatures. Thus, for high temperature
differences between the bodies exchanging energy, radiation becomes the domi-
nant heat transfer mode. This strong nonlinear dependence of the emitted radia-
tion on absolute temperature T as T4 makes even simple conduction problems
extremely difficult to solve analytically.
In the next section, numerical treatment of radiation boundary condition is
discussed with the aid of an example problem.
Example problem Consider an infinite slab of thickness L (Fig. 4.40). Initially,
the body is at a uniform temperature, Ti. Suddenly, one side is exposed to a very
hot environment at temperature Too, the other side being insulated. We wish to
find the temperature distribution in the slab as a function of time.
Assumptions
1. One dimensional heat transfer is taking place because the slab is infinite
in y and z-directions.
It may be noted thatthe coefficient ofTt in Eq. (4.184) will depend on time. This
does not create much difficulty because this term will have to be recomputed at
each time step before going on to the next time step.

4.28 IRREGUlAR GEOMETRY

We have considered so far surfaces parallel to the x and y axes. For such simple
cases, we can usually arrange for certain grid points to lie on the boundaries.
However, there are cases in practice where the surface is curved or irregular and
the boundary does not fallon regular grid points (Fig. 4.41). Hence, special treat-
ments are needed for handling the grid points in the neighbourhood of an irregular
boundary. Two types of situation may arise: (1) When the temperature of the
irregular boundary is specified. (2) When heat flux or normal derivative is
specified on the irregular boundary.
 Numerical Methods for Conduction Heat Transfer 111

superior to Gaussian elimination method. The line-by-line method for faster

convergence is mentioned. The formulation of a three-dimensional problem in
Cartesian coordinates is briefly shown. Subsequently, the finite-difference
formulation of a one-dimensional transient heat conduction problem in a plane
slab is taken up. Three methods of solution, namely, Euler, Crank-Nicholson and
pure implicit are introduced and compared with each other in respect of their
stability and order of accuracy. A physical explanation of the three methods is
also given. The concept of consistency with reference to transient problems is
broached to the readers. For two-dimensional transient problems, the superiority
of Alternating Direction I1J1plicit(AD!) method is then established followed by
an example problem formulation by AD!. An example of a grid independence test
is then given. False transient approach for steady-state problems and ADI method
for three-dimensional transient problems are then discussed. The treatments of
Fig.4.44 Irregular boundary with normal derivative specified on it C(} > 45°) condition at the centre for axisymmetric and non-axisymmetric problems in
cylindrical geometry are given in detail, followed by the mention of the
2. The equations developed for grid points in the vicinity of the irregular application of the same for spherical geometry. Handling of the interface of a
boundaries replace the boundary equations used for rectangular domain. composite body is then detailed. Finally, the treatments of nonlinearities arising
3. If it is desired to obtain T at the boundary grid points, they maytbe out of nonlinear governing differential equation and/or boundary conditions are
computed using linear interpolation between the boundary point desired discussed followed by handling of irregular geometry with the aid of example
and adjacent interior grid points whose values were calculated earlier.
problems.
Conclusion The above discussion on treatment of irregular geometry reveals
that although finite-difference is capable of handling irregular or curved bound-
ary, the procedure is rather tedious as lots of account keeping is necessary for grid REFERENCES
points near the boundary. A better alternative is to use grid generation technique
or finite element method which are better suited for handling irregular geometry. 1. Arpaci, V S, Conduction Heat Transfer, Addison-Wesley, 1966.
2. Gebhart Benjamin, Heat Conduction and Mass Diffusion, McGraw-Hill Inc.,
New York, 1993.
3. Brian, P L T, "A Finite-Difference Method of High-Order Accuracy for the
SUMMARY
Solution of Three-Dimensional Transient Heat Conduction Problems", A.I.Ch. E.
Journal, Vol. 7, 1961, pp. 367-370.
This chapter begins with the finite-difference formulation of the problem of a
4. Carnahan, Brice, H A Luther, and James 0, Wilkes, Applied Numerical Methods,
one-dimensional steady-state heat conduction in a cooling fin. The treatment of
John Wiley & Sons, New York, 1969.
insulated and convective boundary conditions by the use of image-point technique 5. Constantinides, Alkis, Applied Numerical Methods with Personal Computers,
is discussed. Alternative to the image-point technique such as use of polynomial McGraw-Hill Inc., New York, 1987.
fitting and its application is also discussed. The concept of tridiagonal matrix is ( 6. Frankel, S P, "Convergence Rates of Iterative Treatments of Partial Differential
brought forth. The resulting system of equations is solved by Gaussian Equations", Math. Tables Aids Comput., Vol. 4, 1950, pp. 65-75.
Elimination, Thomas algorithm (tridiagonal matrix algorithm or TDMA) and 7. Jaluria, Yogesh, Computer Methods for Engineering, Allyn and Bacon, Inc.,
Gauss-Seidel iterative methods and the superiority of Thomas algorithm for Boston, 1988.
solving systems having tridiagonal coefficient matrix is demonstrated. The 8. Jaluria, Yogesh and Kenneth E. Torrance, Computational Heat Transfer,
concepts of over-relaxation, under-relaxation and optimum relaxation factor are Hemisphere, Washington, D.C., 1986.
introduced with respect to Gauss-Seidel iterative method. Then the formulation of 9. James, M L, G M Smith, and J C Wolford, Applied Numerical Methods.for Digital
Computation with FORTRAN, Scraton, International Text Book Co., 1967.
a two-dimensional steady-state problem with heat generation in a square slab is
10. Myers, Glen E., Analytical Methods in Conduction Heat Tran~fer, McGraw-Hill,
shown. Handling of the comer points is demonstrated. It is also shown that for
New York. 1971.
large system of equations with banded coefficient matrix, G-S iterative method is 11. Patankar, S V, Numerical Heat Transfer and Fluid Flow, Hemisphere,
Washington, D.C., 1980.
5.1 INTRODUCTION
This chapter deals with the calculation of flow field of incompressible (constant
density) fluids by finite-difference methods. To obtain the velocity and pressure
distribution in a fluid flow one has to solve the equations of conservation of
momentum (Navier-Stokes equations) and conservation of mass (equation of
continuity) simultaneously. The important practical applications of numerical
methods for the solution of Navier-Stokes equations are flow around aerofoils,
wings, fuselages and ship hulls, cars, trains and so forth which ultimately results
in calculation of drag over the said objects. Other applications include flows in
pipes and channels, and bearing lubrication. It may be noted that flows may be
associated with heat transfer as in convection heat transfer to be separately dealt
with in Chapter 6.
Three classes of flow problems will be discussed here: (i) Creeping flow (the
limiting case of very large viscosity, that is, very small Reynolds number), (ii)
Boundary layer flow (the limiting case of very small viscous forces, that is, very
large Reynolds number), and (iii) Inviscid flow or frictionless flow (ideal fluid,
J.l = 0). In all three cases, the flow geometry is taken as rectangular. Furthermore,
flow is assumed as laminar and isothermal and viscosity is not a function of
temperature.
It may also be noted that gases may be treated as incompressible fluids (when
the speed is low, Mach number less than about 0.3, and the external heating or
cooling is small). Hence, the study of incompressible viscous flow is not so lim-
ited as it might appear and its applications range from oceanography to lubrica-
tion to aerodynamics (Hughes, 1979).
 Numerical Methods for Incompressible Fluid Flow 121

The pressure field is indirectly linked with the continuity equation. When the
correct pressure field is plugged into the momentum equations the resulting
velocity field satisfies the continuity equation.
Therefore, it should now be clear that calculation of pressure field is the main
problem facing the numerical analvst in solving fluid flow problems. Then
following sections deal with the algorithms for obtaining pressure and velocity
field by two different methods, namely, Stream Function- Vorticity and Primitive
Variables methods.

5.4 STREAM FUNCTION-VORTICI1Y METIIOD

In this method, the difficulty associated with the computation of pressure is

circumvented by eliminating the pressure gradient terms from the momentum
equations by cross-differentiation which leads to a vorticity-transport equation.
This, when coupled with the definition of stream function for steady two-
dimensional situations, is the basis of the well-known stream function-vorticity
method.
Advantages
1. The pressure makes no appearance.
2. Instead of having to deal with the continuity and two momentum equations,
Equations (S.Ia and b) are the Navier-Stokes equations and Eq. (5.2) is the
we need to solve only two equations to obtain stream function and
equation of continuity. For the derivations of the above equations readers are
vorticity.
referred to Schlichting (1979).
Disadvantages
It should be noted that pi denotes the difference between the total pressure and
I
1. Calculation of pressure.
hydrostatic pressure (pressure at rest). This causes the body forces to cancel, as
2. Difficulty in specification of vorticity at a wall.
they are in equilibrium with the hydrostatic pressure.
3. The method is valid for two-dimensional problems as the definition of
5.3 DIFFICULTIES IN SOLVING NAVIER-STOKES stream function applies to two-dimensional flow field. However, the
method can be extended to three-dimensional problems at the expense of
EQUATIONS
considerable amount of complexity using six dependent variables, namely,
1.NonlineQrity A close look at Eqs (S.Ia and b) reveals that the convection part the three components of the vorticity vector and three components of the
of the momentum equations involve nonlinear terms. For example, in Eq. (S.Ia), velocity-potential vector.
the convection coefficientpu is a function of the dependent variableu. However,
this can be treated like the conductivity, k being a function of temperature, T as
discussed in Chapter 4. Starting with a guessed velocity field, one could itera-
tively solve the momentum equation to arrive at the converged solution for the
velocity components.
Therefore, nonlinearity poses no problems as such. It only makes the computa-
tions more involved.
II. Pressure gradient The main hurdle to pvercome in the calculation of velocity
field is the unknown pressure field. The pressure gradient behaves like a source
term for a momentum eq1.tation.But, there is no equation for obtaining pressure.
For a specified pressure field there is no particular difficulty in solving the
momentum equations. So, the challenging task is to determine the correct pressure
distribution.
5.7 INVISCID FLOW (STEADY)

For steady inviscid flow.(u = 0), vorticity is constant along a streamline. There-
fore, vorticity distribution at the inlet calculated from the inlet velocity distribu-
tion can be applied everywhere in the computational domain. Hence, the equation
to solve is Eq. (5.11) for irrotational flow and Eq. (5.12) for rotational flow.
~ can be computed from Eq. (5.3). Once the Laplace equation for stream function
is solved, the velocity distribution in the domain can be calculated from Eqs (5.6)
and (5.7).
For specifying boundary conditions it should be kept in mind that for inviscid
flow, since the viscous terms are neglected, the order of the governing equations
of motion drops from two to one. As a result, only one boundary condition, with
respect to the velocity field, can be satisfied at the boundaries. Thus; slip is al-
lowed parallel to the walls (in contrast with the no-slip condition for viscous
flows) and the normal velocity component is taken as zero. Therefore, the value
of stream function is constant along a wall.
5.8.1 Boundary Conditions for Pressure
It is obvious that to solve Poisson equation for pressure, four boundary condi-
tions, two in each direction are needed. These boundary conditions may be ob-
tained from x- and y-momentum equations and these are generally in terms of
gradients.

5.8.2 Solving for Pressure Distribution in

the Domain and on the Wall
For a steady flow problem, to obtain the pressure distribution in the computa-
tional domain, the Poisson equation for pressure (Eq. (5.25» is solved only once,
that is, after the steady state values of; and Vlhave been computed. For unsteady
flow problems, VI-distributions at given time intervals are used.
5.12.5 A Remedy: The Staggered Grid
The aforementioned difficulties can be circumvented by recognising that all the
variables need not be computed for the same grid points. It is possible to employ
a different grid for each dependent variable.
In the staggered grid, the velocity components are calculated for the points that
are located on the faces of the control volume (shaded) while the pressure is cal-
culated for the regular grid points. For uniformly spaced grid points, the control
volume faces are situated exactly at the· midway between the grid points.
Figure 5.4 s~ows how the locations for the velocity components u and v are to be
defined.
Advantages
1. For a typical control volume shown (shaded with hatched lines), the
discretized continuity equation would contain the differences of the
adjacent velocity components and hence a wavy velocity field would not
be satisfied.
2. The pressure difference between the successive grid points now becomes
the natural driving force for the velocity component located between the
grid points. Hence, a non-uniform or a wavy pressure field will not be
treated as a uniform pressure field and cannot arise as possible solutions.
136 Computer Simulation of Flow and Heat Transfer

corrections and velocity corrections at the neighbours o!lurroJUlding grid

points. These neighbours, would in turn, bring their nli,,.,,.,,, and so on.
Ultimately, the velocity correction formula would involvI thlopressure
correction at all grid points in the calculation domain and thl resulting
pressure-correction equation becomes unmanageable. The omu.Jionof any
term would, of course, be unacceptable if it meant that the ultimate solution
would not be true solution of the discretized form of the momentum' and
continuity equation. It so happens that the converged solution given by SIMPLE
does not contain any error resulting from this omission. In the converged
solution we acquire a pressure field such that the corresponding velocity field
do satisfy the continuity equation.
The words, semi-implicit in the name SIMPLE have been used to acknowl-
edge this omission. This term represents an indirect or implicit influence of the
pressure correction on velocity; pressure corrections at nearby locations can
alter the neighbouring velocities and thus cause a velocity correction at the
point under consideration. We do not include this influence and thus work with
a scheme that is only partially, and not totally implicit.
Since the velocity corrections can be assumed to be zero at the previous itera-
tion step, Eq. (5.44) and Eq. (5.45) can be written as
5.13 COMPUTATION OF BOUNDARY lAYER FLOW

5.13.1 The Concept of Boundary Layer

When fluids of small viscosity such as air or water move rapidly over a solid body
(thus making Reynolds number of the flow very high), the friction between the
fluid and the solid surface causes the fluid to be retarded within a thin region
immediately adjacent to the solid surface. This thin region, in which large veloc-
ity gradients exist, is called boundary layer (introduced by L. Prandtl in 1904).
The region outside the boundary layer, where the forces due to friction are small
and may be neglected, and where the real or perfect fluid approximation is valid
is called potential or inviscid flow region (Fig. 5.9).
 Numerical Methods for Incompressible Fluid Flow 149

variables approach, the concept of staggered grid has been introduced and the
SIMPLE algorithm of Patankar and Spalding (1972) has been described in detail.
The method of solution of biharmonic equation in If! arising in creeping flow
problems has been detailed. Finally, the method of computation of boundary layer
flow over a flat plate and a wedge by similarity method using shooting technique
is discussed. The chapter ends with a treatment of finite-difference methods using
x-y coordinates and also X-OJ coordinates (Von Mises transformation) for flat
plate boundary layer problems.

REFERENCES

1. Anderson D A, J C Tannehill, and R H Pletcher, Computational Fluid Mechanics

and Heat Transfer, Hemisphere, Washington, D.C., 1984.
2. Hughes, William, An Introduction to Viscous Flow, Hemisphere, Washington,
D.C., 1979.
3. Jaluria, Yogesh and Kenneth. E. Torrance, Computational Heat Transfer,
Hemisphere, Washington,D.C., 1986.
4. Patankar, S V, Numerical Heat Transfer and Fluid Flow, Hemisphere,
Washington, D.C., 1980.
5. Patankar, S V, and D B Spalding, "A Calculation Procedure for Heat, Mass and
Momentum Transfer in Three Dimensional Parabolic Flows", Int. J. Heat and
Mass Transfer, Vol. 15, 1972, pp. 1787-1805.
6. Schlichting, H, Boundary-layer Theory, McGraw-HillBook Company, 1979.
7. Shih, J M, Numerical Heat Transfer, Hemisphere, Washington, D.C., 1984.
8. Tadmor, Zehev and Costas G. Gogos, Principles of Polymer Processing, John
Wiley & Sons, New York, 1979.
9. Welch, J E., J H Harlow, J P Shannon, and B J Daly, "The MAC Method: A
ComputingTechnique for Solving Viscous, IncompressibleTransient Fluid Flow
Problems Involving Free Surfaces", Los Alamos Scientific Lab. Rept. LA-3425,
1965.
From Eq. (5.78), it is readily seen that the transverse velocity v does not '
appear. This is a major advantage of Eq. (5.78) over Eq. (5.69) since the
computation of one less unknown will be required in the present case. EXERCISE PROBLEMS
For more details, readers are referred to Shih (1984).
5.1 Consider the inviscid flow between two parallel plates as shown in
Fig. P. 5.1. The fluid enters the channel at a uniform velocity of2m1s. The
SUMMARY velocity at the inflow, as well as that at the outflow, is 5 mls. The inlet and
outlet (at the bottom plate) are each 4 cm wide, and their centrelines are
This chapter deals with the method of computation of flow field of incompressible 12 cm apart. If the vertical distance between the plates is 20 cm and the
fluid by finite-difference techniques. Two types of flow problems have been exit channel is 1 cm high, obtain the streamlines and the velocity distribu-
discussed namely.viscous flow (including creeping flow and boundary layer flow) tions (u and v) for this flow, assuming it to be two-dimensional. The
and inviscid flow, with greater emphasis placed on the former. Both stream velocity at the exit can be taken as uniform assuming it to be a long narrow
function-vorticity and primitive variables approaches have been discussed in passage.
detail with the mention of advantages and disadvantages of each. In the primitive
152 Computer Simulation of Flow and Heat Transfer

5.5 Solve the problem 5.4 by stream function-vorticity approach for air
instead of oil flowing through the enclosure. Note that creeping flow
approximation will no longer be valid here. Why?

6.1 INTRODUCTION
The convection heat transfer process can be divided into two broad categories
(a) when the flow is generated due to an external forcing agent, such as fan,
pump, natural breeze, or the motion of a heated object, the process is called forced
convection. This is particularly important circumstance in industry, since the flow
is induced by means of an external agency in most thermal systems, such as heat
exchangers; (b) when the fluid motion arises simply because of density differ-
ences, caused by temperature differences in a body force field, such as gravity,
the process is termed natural or free convection. A heated body cooling in still air,
the free rise of heated buoyant air due to a fire, circulation of water in a lake,
dissipation of heat from the coil of a refrigeration unit to the surrounding air, are
some of the examples of natural convection.
There is also another class of convection known as mixed convection in which
both natural and forced convection effects are significant and neither mechanisms
may be neglected. Manufacturing processes such as wire drawing and extrusion
belong to the aforementioned category.
There is an important difference between forced and natural convection. In the
former case the flow is externally imposed and is often independent of the
temperature field. The flow field can thus be obtained independent of the heat
transfer processes and then used in the determination of the temperature field. In
natural convection, on the other hand, the flow and heat transfer are coupled,
since the flow itself arises due to the temperature difference in the fluid. As a
result, the flow field cannot be obtained independent of the temperature field and
the two must be considered simultaneously. This complicates the problem
considerabJy, and thus even the simplest problem usually demands numerical
.Wlb.ods of solution. In forced convection too, if fluid properties are functions of
1S4 Computer Simulation of Flow and Heat Transfer

the flow field is not independent of temperature field and the two are linked,
giving rise to a complexity similar to that in free convection.
In the following sections, we shall discuss numerical methods for forced
convection in which fluid properties are assumed to be constant, i.e. the tempera-
ture variations within the fluid are not large enough to affect its properties. First,
we will consider the class of problems commonly ~wn as convection-diffusion
type of problems. Next, the thermal boundary layer flows will be dealt with.

6.2 CONVECTION-DIFFUSION (STEADY,

ONE-DIMENSIONAL)

Let us now examine the simultaneous effects of steady axial convection and diffu-
sion in slug-flow with velocity, Uo, subjected to specified temperature To for
x $ 0 and Tl for x ~ L (Fig. 6.1). We seek the steady state temperature T(x) by
numerical means.
lS8 Computer Simulation of Flow and Heat Transfer

Now, since ~ee = Pe &X, &x has to made smallet because Pe is large. In other
words, grid is to be refined.
But, refinement of grid is not computationally efficient, as it will require more
computer time as a large number of equations will then have to be solved.

6.2.4 Exponential, Hybrid and Power-law

Schemes
There are other schemes such as exponential schemes where the discretization
equation is obtained directly from the exact solution. This is to avoid the defects
of the upwind·scheme enumerated in the previous article. The method is attractive
in terms of accuracy, but the computation of exponentials is time consuming. The
hybrid scheme combines upwind and central differencing of the convective terms
to achieve the stability of the upwind method and some of the better formal accu-
racy of the central differencing. The scheme of discretization of the convective
term is changed from central to upwind when Pee;;::2. A much better accuracy is
obtained by the use of a piecewise power-law scheme, discussed by Patankar
(1980). A power-law scheme is employed in the range of 0 ~ Pee ~ 10 so that
central differencing results at Pee = 0 and upwind differencing for Pee> 10. The
power law scheme is a computationally efficient approximation of the exponen-
tial method.
However, the power-law scheme is strictly valid for the steady one dimen-
sional case and cannot be extended to two or three dimensional or transient prob-
lems. In most practical situations, one has to deal with multi-dimensional convec-
tion-diffusion problems and therefore, a need for a more general method has
arisen. In Chapter 7, a method known as operator-splitting (O-S) algorithm
employed to alleviate the defects of the upwind scheme at high Peclet numbers
will be discussed. The O-S method is essentially an alternative method to upwind
differencing.

6.3 CONVECTION-DIFFUSION (UNSTEADY,

ONE-DIMENSIONAL)
Let us now consider a problem in which a hot fluid with unit velocity is intro-
=
duced suddenly (at t 0) at the inlet of a flow domain (Fig. 6.4). The cold fluid in
the physical domain (0 <x < 00) is initially at a constant temperature. If the Peclet
number is large, a thermal front whose thickness is small will move through the
flow region with unit velocity. As Pe ~ 0, the conduction (or diffusion) limit is
reached and hence the front will be diffuse and the energy transport is due to the
molecular action of thermal diffusivity.
The non-dimensional differential equation is,
 Note that the auxiliary variables Ej' Gj are (2 x 2) matrices and Fj is a (2 x 1)
vector, implying the solution of system of four and two algebraic equations.
To sum up, the main features of the implicit Keller-box method are:

• Only slightly more algebra than in Crank-Nicholson method.

• Second order accuracy with arbitrary (non-uniform) x andy-spacings.
• Allows very rapid x-variations.
• Allows easy programming of the solution of large number of coupled
equations.

6.7 TRANSIENT FREE CONVECTION FROM A

HEATED VERTICAL PLATE

6.7.1 Introduction
Free (or natural) convection currents in the vicinity of a heated object such as a
room heater placed in air are familiar phenomena. The density of the fluid (in this
case air) decreases, resulting in a buoyancy effect causing circulation of the fluid.
Here the velocity and temperature distributions are interrelated; the temperature
distribution, in effect, produces the velocity distribution. Figure 6.12 shows the
 Numerical Methods for Convection J:Ieat Transfer 175

three-dimensional flow, one must use the primitive variables (u, v, T) approach to
arrive at a solution.

SUMMARY

In this chapter, finite-difference formulations for both forced and natural convec-
tion problems have been discussed. The chapter begins with a detailed study of
the possible numerical methods for solving 1-D steady convection-diffusion prob-
lems. It is shown that upwind scheme is superior to central differencing when the
Peelet number of the flow is high. The concept of false diffusion arising in
upwinding is also introduced. Brief mention is made of exponential, hybrid and
power-law schemes. Also, finite-difference formulations using upwinding is
described with respect to unsteady, one-dimensional and two-dimensional con-
vection-diffusion problems. Next, computation of thermal boundary layer flows
(both external and internal) is taken up. Both explicit and implicit methods of
801utionare discussed in detail for computation of temperature field of fluids of
Pr ~ 1 flowing over a flat plate maintained at constant temperature. For computa-
tion of temperature field of liquid metals (0.005 < Pr < 0.05) flowing between
parallel plates-the upper plate insulated and the lower plate at constant tempera-
ture, three methods namely, explicit, Crank-Nicholson and Implicit Keller-box
,e described in detail. The concept of block tridiagonal matrix algorithm is
~roached. Under the topic of natural convection, the problem of the solution meth-
odology of the transient free convection from a heated vertical plate is given in
cJetail using both explicit and implicit schemes. Finally, only the qualitative
aspect of the method of solution of the problem of free convection in enclosures is
mentioned briefly.
!.,.

REFERENCES

1. Arpaci, V S and P S Larsen, Convective Heat Transfer, Prentice-Hall, New Jersey,

1984.
2. Carnahan, Brice, H A Luther, and James 0 Wilkes, Applied Numerical Methods,
John Wiley & Sons, New York, 1969.
3. Hellums, J D and S W Churchill, "Transient and Steady State, Free and Natural
Convection, Numerical Solutions: Part 1, The Isothermal, Vertical Plate",
A.I.Ch.E. Journal, Vol. 8, 1962, pp 690-692.
4. Jaluria, Yogesh and Kenneth E Torrance, Computational Heat Transfer,
"., Hemisphere, Washington, DC, 1986.
,0 5. Muralidhar, K and P S Ghoshdastidar, Computer Methods in Heat Transfer, ISTE
(Indian Society for Technical Education) Course Package, New Delhi, 1992.
6. Patankar, S V, Numerical Heat Transfer and Fluid Flow, Hemisphere
Washington, DC, 1980.
'7~1 INTRODUCTION

In this chapter the readers are introduced to two new methods, namely (i) applica-
tion of ADI method to solve 2-D transient heat transfer from a rectangular, sym-
metrically trilaminated composite fin, and (ii) Operator-splitting algorithm which
is an alternative to popular upwind scheme to solve convection-diffusion prob-
lems. First, we shall discuss method-(i) and then, method-(ii).

'1.2 NEW MEmOD (I): APPLiCATION OF ADI

MEmOD TO SOLVE mE PROBLEM OF 2-D
TRANSIENT HEAT TRANSFER FROM A
STRAIGHT COMPOSITE FIN

7.2.1 Application of Composite Fin

It is well known that a fin is used to augment heat transfer between an object and
its surrounding fluid. Under certain conditions, the fin is in contact with high
temperature or corrosive fluid. A layer of high strength or corrosion resistant
J;Ilaterialis coated on both sides of the fin to withstand harsh environment (Ju et
fll. 1991).

'.2.2 Background
tit Chapter 4, the Alternating Direction Implicit (ADI) scheme (which is a direct
or non-iterative numerical method) is discussed with respect to the solution of
transient multi-dimensional conduction problem in bodies having isotropic and
liomogeneous properties. The problem of applying ADI method to composite
bodies is the interface. Ghoshdastidar and Mukhopadhyay (1989) were the first
to successfully apply ADI method to obtain solutions for transient heat conduction
in a rectangular composite fin thus overcoming the problem posed by the
7.3.4 Results
Figures 7.4 and 7.5 show the plots of temperature distribution as a function of
distance for various time levels and two Peclet numbers Pe = 0.1 aild 100. Typi-
cal values of Llx (=I1t) used in these calculations are 0.1 and 0.05. The as results
have been compared with analytical, central difference and upwind results
(Muralidhar and Ghoshdastidar, 1992).
The agreement between the numerical predictions and the analytical result is
very good atPe = 0.1. AtPe = 100, only the as results matches with the analyti-
cal solution closely. The upwind method fails due to false diffusion (see Chapter
6) inherent in its formulation. The central difference method fails due to matrix ill
conditioning and loss of diagonal dominance at high Peclet numbers. The as
algorithm is uniformly accurate at low as well as high Peclet numbers.

7.3.5 Concluding Remarks

The operator-splitting algorithm has been further satisfactorily tested f~r
transient one-dimensional problems including radioactive decay, energy sources
and a non-homogeneous initial condition, and for model two-dimensional
 Special Topics 1 New Methods 189

problems. The algorithm has also been successfully applied to solve engineering
problems such as Adsorption-Desorption problems in porous region, numerical
modelling of enhanced oil recovery using water-injection method and flow past a
discrete protruding heater in a vertical channel (Ghoshdastidar and Muralidhar,
1994).
In conclusion, it can be said that OS algorithm is a viable alternative to upwind
scheme so far as the convection-diffusion problems in heat transfer are concerned
and can be applied effectively to solve practical problems.

SUMMARY

This chapter introduces the readers to two new methods of numerical solution
(developed in recent years) of important heat transfer problems of practical
significance in some detail. They are: (i) application of ADI scheme to the
solution of transient heat transfer problem in a straight composite fin, and
(ii) alternative to upwind scheme-Operator-Splitting algorithm to solve
convection-diffusion problems.

REFERENCES

1. Ghoshdastidar, P S and A Mukhopadhyay, "Transient Heat Transfer from a

Straight Composite Fin: A Numerical Solution by ADI", International
Communications in Heat and Mass Transfer, Vol. 16, No.2, 1989, pp. 257-265.
2. Ghoshdastidar, P Sand K Muralidhar, "Numerical Solution of Convection-
Diffusion Problems without Using Upwind Method", Project Report of DST
(Department of Science and Technology) Project No. III. 5 (l41)/88-ET, New
Delhi, 1994.
3. Ju, Yi-Hsu, Yen-Ping Shih, and Chen-Yaw Chin, "Design of a Trilaminated
Rectangular Fin", International Journal of Heat and Mass Transfer, Vol. 34,
No. 4/5, 1991, pp. 1097-1104.
4. Muralidhar, K and P S Ghoshdastidar, Computer Methods in Heat Transfer, ISTE
(Indian Society for Technical Education) Course Package, New Delhi, 1992.
 Computer Modelling of Some Important Industrial Problems Showing
Applications of Radiation in Enclosure, Combined Modes of
Heat Transfer, Non-Newtonian Flow and Phase Change

8.1 INTRODUCTION
This chapter presents applications of various numerical methods to model some
important industrial problems. In each case, considerable details have been given
to make the readers appreciate the modelling techniques and interpretation of
results. For some cases, a separate nomenclature section has been added at the
end of the discussion.

8.2 TRANSIENT COMBINED MIXED CONVECTION

AND RADIATION FROM A VERTICAL
ALUMINIUM FIN (GHOSHDASTIDAR
AND RAJU, 1992)
Introductory Remarks This discussion presents a numerical solution of
transient conjugate heat transfer from an isolated vertical rectangular air-cooled
aluminium fin. Heat transfer by both mixed convection and radiation is
considered. Mixed convection effect should be appreciable for low speed air flow
over the fin, e.g. when the fin is stationary and wind is blowing over it at a low to
moderate speed. Radiation heat transfer mode is important for large temperature
difference between the fin and the surrounding as well as for high emissivity fin
material and for low speed air-flow. Consideration of transients is essential as
during start-up and shut-down the heat transfer rates vary significantly. Conjugate
heat transfer should be taken into account as it gives rise to non-uniform heat
transfer coefficient the consideration of which is more realistic in contrast with
the assumption of uniform heat transfer coefficient in conventional fin theory.
The objectives of the present study are (i) to develop a good numerical
procedure to handle transient conjugate heat transfer from a vertical fin when
radiation is taken into account and (ii) to see how important is the radiation mode
 Special Topics II 193

The time-step, Ll't', the grid size !::J( and LlY have been taken as 0.02, 0.1 and
0.002 respectively. Maximum values of X and Yare 1.0 and 0.03 respectively.
Uniform grid spacings have been used in X and Y directions. The grid sizes and
time-step are chosen after having performed sufficient numerical
experimentations. An 11 x 16 grid system has been used.
Algorithm Due to non-linearities and mutual coupling, the governing equations
are solved simultaneously and iteratively. The iterative procedure involves an
overall iteration loop together with sub-iteration loops.
The mixed convection loop (see Fig. 8.2) computes convection heat transfer
coefficient, 'he' as a function of position and time which means a set of values of
'he' as a function of temperature has been obtained. This is required to solve the
fin heat conduction Eq. (8.4b) which includes radiation effect also.
To solve Eq. (8.4b), the future fin temperature distribution is assumed first.
This is required to calculate the radiation heat transfer coefficient, h,. To make an
educated guess, the future fin temperature distribution is assumed to be the one at
the end of the time-step ('t' + Ll't') for the case of only pure mixed convection.
Equation (8.4a) and Eq. (8.4b) are solved by Tri-diagonal Matrix Algorithm
(TDMA).
The algorithm, in the form of a flow chart is shown in Fig. 8.2.
The numerical results have been obtained using the data listed in Table 8.1.
The fin emissivity value at the base temperature Tb is used. all other properties
are evaluated at the average of the base and the ambient temperature, i.e. at
(Tb + T,,J/2. The computations have been carried out on an IBM PC-XT compat-
ible. Average execution time per computer run (i.e. for a particular set of param-
eters) is approximately eight minutes.

Table 8.1 Input Data to the Program

Fin (Aluminium)
Length of the fin (L) =0.1 m
Base temperature (Tb) = 773 K
Thermal conductivity (ks) = 204 W/m-K
Emissivity (e) = 0.049
Fluid (Air)
Prandtl number (Pr) = 0.7
Temperature (T~) = 300 K
Reynolds number (Re) = 50,000
Thermal conductivity (kf) = 0.045 W/m-K

Results and Discussion Figure 8.3 shows the temperature distribution of the fin
at different times. It is clear that with the progress of time, the temperature of the
fin increases which is also expected from the physical point of view. Numerical
results show that aow / ax is not equal to zero at the tip of the fin for larger tim~s.
This error can be removed by making finer grid in the fin.
 Fig. 8.5 Temperature vs. time plots at X = 0.5 of the fin
for CCP = 5 and for different Gel Rfl values

gradient is steeper indicating more heat transfer from the fin. This is expected as
fin temperature as observed from the plots. The dashed lines are the results of the increased Gr/Re2 implies increased vertical velocity of the fluid due to larger
conventional fin theory (Sunden, 1983). The results of the conventional fin theory buoyancy force.
and the present numerical method show similar trend, although the former gives Figure 8.11 indicates vertical velocity profiles of the fluid in the steady state
higher temperature in the fin. It may be noted here that while heat transfer coeffi- for Gr/Re2 = 0, 2 and 10. The velocities are higher for greater Gr/Re2 as
cient is constant over the entire fin surface in the conventional fin theory the same increased Gr/Re2 means buoyancy effect is more. For Gr/Re2 = 10, the velocity
is non-uniform in the present study. =
profile reaches a peak (U 1.15) and then decreases gradually to U 1.0, the=
In Fig. 8.8, steady state fin temperature profile for pure forced convection reason being the flow is aided much more by the buoyancy force in comparison
=
i.e. Gr/Re2 0 at different CCPs are presented. The fin temperatures are higher with the case of Gr/Re2 = 2.
than those for ~e cases shown in Fig. 8.7. This is because when Gr/Re2 is Figure 8.12 provides distribution of the local heat flux for CCP = 5 and
positive, the buoyancy force assisted by the main forced flow increases the mixed Gr/Re2 = 2. At the leading edge the heat flux is still infinite (since X = 0 is a
convection heat transfer resulting in lower fin temperature. singular point). With increasing X, the heat flux decreases at first, reaches a mini-
Figure 8.9 shows the transient fluid te~peratures in the y-direction (atX = 0.5) mum and then increases to its value at X = 1.0.
at different times. It is seen that at greater times .the temperature gradient at the Table 8.2 shows that the contribution due to radiation is much smaller com-
wall is steeper indicating larger heat transfer. That is, initially heat transfer rate is pared to that due to mixed convection. This is because, the temperature difference
less and then gradually it increases. in this case between the fin and the fluid is not extremely high, the emissivity of
=
Figure 8.10 shows the steady state fluid temperature for Gr/Re2 0 and 2 at the fin material is low and also Reynolds number is high (Re = 50,000). The
=
X 1.0. With higher Gr/Re2, the fluid temperature is less and the temperature radiation effect can be significant for very high temperature difference, low speed
8.3 HEAT TRANSFER IN ROTARY KILN
REACTORS

8.3.1 Introductory Remarks

An industrial rotary kiln consists of an open ended long, refractory lined cylindri-
cal shell mounted at aslight incline from the horizontal plane. The solid raw
material ("raw mixture") is fed into the shell at the upper end and during the
process becomes transferred to its lower end, where it is discharged.
In the majority of high temperature kilns the required reaction temperature of
the raw material is achieved by gas or oil combustion in the kiln (internally fired
kilns) above the solid, usually by employing a burner at the end of the kiln and
producing a flow of burning gas countercurrent to the solid movement. The rotary
kilns have low relative load (only 10 to 20% of the kiln's cross-section is occu-
pied by solids). Industrial rotary kilns are often very large (up to 220 m in length
206 Computer Simulation of Flow and Heat Transfer Special Topics II 207

and 5 m in diameter), however, those employed in the chemical industry are the calcining reactions (CaC03 ~ CaO) begins (Spang, 1972). The present
smaller (Helmrich and Schtigerl, 1980). study is limited to the preheat zone where no chemical reaction occurs.
The rotary kilns are used mainly in the cement and metallurgical industries. Its Figure. 8.14 (a) shows the nomenclature of the kiln. The angle, ais called the
other applications include production of chemicals and fertilizers such as phos- fill-angle which shows the region of volume containing the solid. Figure 8.14(b)
phates and polyphosphates, carbonates, zinc white, barium sulphide (BaS), soda shows the heat transfer processes occurring in the kiln.
ash and alumina. The rotary kilns are also used for burning of toxic or non-toxic
waste.
Some of the important advantages of rotary kiln reactors as compared to other
types, such a fluidized bed reactors are simple design and high flexibility, co or
countercurrent gas flow, easy adaptability to changing operating conditions by
varying the length, diameter, inclination, speed of rotation, etc. Processes with
highly corrosive reaction mixtures can be carried out in rotary kilns. Very high
temperature processes can be performed in rotary kilns and materials of varying
consistencies can be handled easily. The main disadvantages are low space-time
yield, inefficient energy utilization and difficult control of temperature and com-
position of the reaction mixture in the central zones (Helmrich and Schiigerl,
1980).
In spite of the diverse uses of rotary kilns in industry, very few systematic
theoretical models of the transport processes in the rotary kilns exist in literature.
There is an urgent need for optimization of rotary kilns.
In this section, two applications of rotary kilns are considered, namely (i)
rotary kiln dryer with applications to the non-reacting zone of a cement rotary
kiln and (ii) rotary kiln solid waste incinerator. The first application involves heat
transfer without chemical reaction while the second one involves heat transfer
with chemical reaction.

8.3.2 Heat Transfer in the Non-Reacting Zone

of a Cement Rotary Kiln
Introduction In a rotary kiln dryer the wet solids are fed into the upper end of
the cylinder and during the process, are dried and heated by the countercurrent Fig. 8.14 (a) Nomenclature of the kiln. (b) Schematic diagram showing
flow of the hot gas. Finally, the solid is transferred to the lower end where it heat transfer processes in a rotary kiln
reaches the desired temperature and is discharged.
In the present study, the kiln can be divided into three sections. In the first PROBLEM FORMULATION
section, the wet solids are heated to the boiling point of the entrained liquid. In the
second section, the liquid evaporates at constant temperature until the feed is Thermal Radiation Among Hot Gas, Refractory Wall and Solid
completely dry. In the third section, the solids are heated to some specified tem- Surface Heat is exchanged among the hot gas, the inner wall of the kiln and the
perature and then are discharged from the kiln. solid surface by radiation as the gas temperature is very high. The wall is divided
For modelling, the non-reacting zone of a cement kiln is considered. The raw into surface elements as shown in Fig. 8.15. Each axial segment (of size equal to
materials fed into the kiln contain calcium carbonate (CaC03), silica (Si02), 1 m) of the refractory surface is divided into fifteen surface elements of equal
shale (AI203) and iron ore (Fe203). These ar~ ground to a very fine powder and size. The solid surface is divided into five surface elements. The temperature of
mixed according to the type of cement being made. Upon heating by hot gases, . the solid surface element and the gas element are assumed to be uniform in each
various reactions occur. There are three important zones: the preheat zone, the ..axial segment.
calcining zone, and the burning zone. In the preheat dry zone, the free water is ~. The surfaces are assumed diffuse and gray. It is assumed that the surface ele-
evaporated and the solid material is heated to,the point (around 1155 K) where .~ ments exchange radiation only with surfaces of the same axial seQment
.t
 Special Topics II 209

Energy Transport in the Solid Mass and energy balances on an element of

solid contained in an axial segment either in the first or in the third section of the
kiln gives the expression for Ts, z +.1z while the same performed on an element of
the wet solids contained in an axial segment in the second section of the kiln gives
the expression for the rate of evaporation, mv' It may be noted that the tempera-
ture of the solids in the second section remains constant at the saturation tempera-
ture of water. It is assumed that water is always available on the solid surface.
The end of the second section is indicated where the cumulative mv is equal to the
total predetermined amount of water to be evaporated per second.
Energy Balance in the Gas Mass and energy balances on an element of the
hot gas contained in an axial segment in the first or third section of the kiln gives
the expression for Tg, z+62 while the same performed on an element of the hot gas
contained in an axial segment of the second section of the kiln gives another
expression for Tg, z+62' In the second section, the superheating of the vapour is
taken into account.

METHOD OF SOLUTION

A computer program is developed to obtain numerical results for the present prob-
lem. Input data are shown in Table 8.3. Finite difference techniques are used and
steady-state thermal conditions are assumed. A grid independence test has been
done. False Transient approach is used to solve the wall conduction equation. The
solution is initiated at the inlet of the kiln and proceeds to the exit. In the second
section of the kiln, the change in shape factors due to removal of water from the
solids is accounted for. The output data consist of refractory wall temperature,
solids temperature, gas temperature, the individual lengths of the first, second and
third section of the kiln and the total kiln length.

RESULTS AND DISCUSSION

Validation with Earlier Experimental and Numerical Results It is

apparent from Fig. 8.16 that the simulation has accurately predicted the gas and
solid temperature profile in the kiln, the fraction of the kiln used for water
evaporation and the fraction used for heating the solids.
However, the kiln length predicted by the present model turns out to be
22 metres as compared to 25.5 metres predicted by Sass (1967), and 27.3 metres
of the actual kiln. The deviation of the present predicted length from that of the
actual kiln may be attributed to uncertainties in some of the input data such as the
composition of the hot gas, contact heat transfer coefficient, the ambient tempera-
. ture outside the kiln, the convective heat transfer coefficient from the kiln to the
surrounding air, the speed of rotation of the kiln, the kiln inclination angle, the
inlet fill-angle and the emissivity of the solid because the aforesaid data are not
specified in Sass (1967) explicitly and therefore, are obtained in the present study
212 Computer Simulation of Flow and Heat Transfer
Special Topics II 213
in the third section of the kiln the axial gas temperature and dT /dz are higher also because of the fact that for the case of higher solid mass flow rate, greater
while axial solid temperature with respect to per cent kiln length is lower. This amount of dry solids have to be heated.
can be explained as follows. For higher amount of water in the feed, the wet solid Figure 8.19 shows that with higher gas flow rate the axial gas temperature is
will have to traverse a greater distance to be completely dry. This means that lower while the reverse is true for the axial solid temperature. This is because the
although same amount of dry solid will have to be heated to the same exit tem- higher gas flow rate implies that the gas residence time is less and hence the kiln
perature and the fill angle remains same in the third section no matter whether the has to be longer. This means the solid has to stay for a greater period of time in the
solid contains more water or less, the total heat transfer to the solid in the last kiln and therefore, solid temperature is high. The gas, on the other hand, loses
section of the kiln for the solid with high water-content will be greater than that more heat to the solids and hence the gas temperature is low.
for solid with low water content, the reason being exposure of the dry solid to
hotter gas in the former case as the heating starts closer to the outlet of the kiln
which is the inlet for the hot gas. Therefore, to achieve the requisite exit tempera-
ture the dry solid will have to travel a smaller length of the kiln and hence dT/dz
is higher while the axial solid temperature at the same per cent kiln length is
lower for the obvious reason.
Figure 8.18 shows that for higher solid mass flow rate, in the third section of
the kiln the axial gas temperature is higher while the axial solid temperature is
less. This can be explained by the fact that for larger solids mass flow rate the fill
angle increases which in turn reduces mean beam length as the gas volume
decreases and hence the emissivity of the gas decreases. This gives rise to lower
loss of heat by the gas by radiation and hence the gas temperature is high. The
solid temperature is less as the heat transfer by radiation to the solid is less and

Fig.8.19 Axial temperature distribution of the solids and the gas for
different mass flow rates of the hot gas

It is interesting to note the very high inlet gas temperature that is required for
the lowest gas flow rate (3 kg/s). Using such hot inlet gas would be uneconomic
although it would drastically reduce the kiln length. It is also found (the graphs
not shown for the sake of brevity) that the axial gas temperatures remain more or
less independent of the angle of inclination (the present range, 2°-5°) of the kiln
or rotational speed (the present range, 1 - 10 rpm).
This is expected as the gas flow is not a function of the aforesaid parameters.
However, the axial solids temperature is greater for higher angle of inclination
of the kiln or for higher rotational speed. Increase in the kiln inclination angle or
the rotational speed increases the predicted kiln length and vice-versa.
 Fig.S.21 Circumferential refractory surface temperature variation

emerges from the solid waste. The temperature increases sharply when the wall is
near the flame surface but rises less rapidly as it gets farther from the flame. The
peak temperature is highest just before being covered by waste. The refractory
surface temperature decreases while covered by waste. The maximum refractory
temperatures are 1820 K and 1970 K for a equal to 180 and 120 degrees, respec-
tively. The maximum temperature is increased when ais greater due to the greater
flame surface area.
Figure 8.22 shows the radial temperature distribution in the refractory. These
temperatures were calculated at the inlet and with a equal to 180 degrees. Tem-
peratures are shown for () equal to 6, 174 and 330 degrees. At () equal to 6 de-
grees, the wall has just emerged from the waste and consequently the surface
temperature is near the minimum. At ()equal to 174 degrees, the maximum sur-
face temperature of 1780 K is reached. At ()equal to 330 degrees, the surface has
been cooled to the minimum by the waste. Except for the thin layer approximately
0.01 m thick, the temperature is seen to remain essentially constant. Conse-
quently, increasing the refractory wall thickness would not effect the thermal
characteristics of the kiln.
 Fig. 8.22 Temperature variation in refractory wall.

The solid waste temperatures are shown in Fig. 8.23 for a equal to 120 and
180 degrees. The temperature increases almost linearly in the entire length of the
kiln. However, the rate of temperature increase is greater for a equal to 120
degrees than when a is equal to 180 degrees. The two curves are terminated at
the point where the mass flow rate of the solid waste is equal to zero. The maxi-
mum temperature of the solid waste for a equal to 120 degree is 560 K and
occurs at 3 m from the kiln inlet. The maximum temperature for a equal to 180
degrees is 610 K and occur at 7.7 m.
Figure 8.24 shows the mass flow rate of the solid waste as a function of axial
position. The mass flow rate remains nearly constant near the inlet since the
volatalization rate is low. A rapid decrease occurs at approximately 2.0 m for the
120 degrees fill angle and at 5 m for the 180 degrees fill angle. The mass flow rate
Fig. 8.24 Mass flow rate of solid waste as function of axial location.
is reduced to zero at 3.0 m and 7.7 m for the 120 and 180 degrees fill angles
respecti vely.
material can be investigated. Finite difference techniques are used to mode the
Summary of Results Heat transfer characteristics of the rotary kiln incinera- heat conduction in the refractory wall and the energy transport by the solid waste
tor was modelled so that the interaction of refractory wall conduction, flame material. The absorption factor technique (see Appendix D for details regarding
radiaton, refractory surface radiation, and chemical kinetics of the solid waste this method) is used to model the thermal radiation heat transfer between the
220 Computer Simulation of Flow and Heat Transfer Special Topics II 221

refractory surface and flame. Numerical results are presented which demonstrate (relatively shallow) of uniform channel depth where the molten material is further
the capabilities of the computer model. The results of primary importance for heated accompanied by further increase in pressure and the material is subjected
design purposes are the refractory wall temperature variation and the length of to high shear, thus enhancing the mixing. However, it may be noted that pressure
kiln required to completely burn the wastes. Each of the parameters listed in rise in the metering section will occur when die opening at the exit of the extruder
Table 8.4 effects incinerator performance and could be investigated. is narrow. For a widely open valve, the pressure will fall.
In the present example, a detailed numerical study of the thermal transport
s.4 MODEllING OF THERMAL TRANSPORT processes in the metering section of a single screw plasticating extruder process-
PROCESSES IN SINGLE-SCREW PIASTICATING ing low density polyethylene (LOPE) has been performed. The computations have
EXTRUDER WITH APPLICATIONS TO been carried out for a uniform barrel temperature distribution and conducting
POLYMER AND FOOD PROCESSING screw. The plastic granules are assumed to be completely melted at the inlet of the
metering section of the screw extruder.
In this section, a detailed numerical modelling of polymer processing in single-
screw extruders is taken up first and then a brief mention is made about various PROBLEM FORMULATION
aspects of food processing including its basic difference with respect to polymer
processing. Based on the computer models, design and optimization of screw ex- Physical Description of the Model The primary difference between a quasi
truders can be achieved. two-dimensional model and a quasi three-dimensional model is the inclusion in
the latter of the cross convection terms, i.e., thermal convections normal to the
S.4.1 A Quasi Three-Dimensional Computer screw flights and the base of the screw channel. Similar to the quasi two-dimen-
Model of Flow and Conjugate Heat Transfer sional analysis, the flow is assumed to be hydrodynamically developed but
in the Metering Section of a Single-Screw thermally undeveloped. Both the x-boundaries are assumed to be insulated
considering that a negligible amount of heat is conducted between the melt and
Plasticating Extruder
the uncooled screw (Fenner, 1977). Also, two-dimensional conduction is assumed
Introduction Screw extrusion is widely used for the manufacture of plastics, to be present in the x-y plane of the screw body with both the boundaries in the
polymers, pharmaceuticals and food products. Typically, a single-screw t-direction acting as insulated considering that the screw body is sufficiently thick
plasticating extruder consists of a screw rotating uniformly inside a barrel main- in that direction. The conjugate heat transfer at the fluid-screw interface is mod-
tained at a fixed temperature (Fig. 8.25). At one end, there is a hopper through dIed by assuming the screwbody behaving like a body in which heat transfer is
which plastic granules are fed in. taking place in a very thin layer below the screw surface while the rest of the
,crew is adiabatic as the temperature drop is expected to be very small in the
,crewbody. As in earlier two models, the non-dimensional channel length is taken
is 250.

W:ajorAssumptions (a) For ease of analysis, the coordinate system is fixed to

:he screw (Fig. 8.26) and thus the barrel moves in a direction opposite to the
;crew rotation.
(b) As the screw channel is shallow, the curvature effects have been neglected,
md therefore the screw is treated as unwound (see Fig. 8.26).
(c) The fluid is considered to be homogeneous and flow and heat transfer
:teady and quasi three-dimensional.
(d) As convective effects are small in comparison to viscous effects because of
ugh melt viscosity and small hydraulic diameter (implying low Reynolds number
low) creeping flow approximations are made for the conservation of momentum
Schlichting, 1979).
(e) The leakage across the screw flights have been neglected.
 Special Topics II 225

Results and Discussions In this section, the results based on the quasi three-
dimensional model for the screw configuration and input data of Karwe and
Jaluria (1990) are presented. If is important to note that the screw surface tem-
perature in Fig. 8.29 and Fig. 8.31 is the integrated average over x at the screw
surface.
The input data used are that for a typical non-Newtonian fluid (low density
polyethylene) following power-law behaviour having n = 0.5. The results have
been obtained for G = 10.0, Pe = 5000.0, /3 = 1.61, qv = 0.3, ¢J = 16.54°, kf= 0.3
W/m. K and k, (for steel) = 45 W/m-K. The value of ks has been taken from
Holman (1981). A typical WI H equal to 5 (Palit and Fenner, 1972) has been used
in this study. A non-uniform grid gradually increasing in the / -direction away
from the walls has been used. In the y *-direction, i.e. along the height of the chan-
nel, grid spacing is uniform. The computations have been carried out over
23 x 17 grids and 23 x 21 grids in the x-y plane (i.e., the cross-section) of the
channel and the screw body respectively. HJH is taken as 10. ~y;and ~Y*s are
0.0667 and 0.5263 respectively. L1z*is 2.5 and the computations stop at
z* = 250.0. Sufficient numerical experimentations have been performed for the
choice of optimum number of grid points .
.Velocity Vectory Plots and Isotherms in the x-y Planes of the
Channel Figure 8.27 shows the velocity vector plots in the x-y plane at four
downstream locations. The plots clearly show the circulatory nature of the flow in
the cross-sectional plane ofthe extruder. This circulation produces the cross ther-
mal convection enhancing the overall mixing process. It is also noticed that v *
velocities are negligible except near the right and left boundary walls meaning a
very weak flow and thermal convection normal to the channel base.
Figure 8.28 shows isotherm plots in thex-y plane at four locations in the down-
~channel direction. At z * = 12.5, i.e. near the inlet of the channel, it is seen that a
~Jarge core of relatively cold fluid is sandwiched between hot fluids in the upper
and lower region of the channel cross-section. Near the right and left boundaries,
temperatures are large because the sides are assumed to be insulated.
At higher Z*,i.e. atz* = 82.5, a dramatic rise in the fluid temperature is noticed.
This is due to the fact that more heating of the fluid is now taking place as the
fluid is exposed to larger surface of the heated barrel and the effect of viscous
heat generation is now more pronounced. The temperature gradient in the y-direc-
tion as well as in the x-direction decreases which means that compared to that at
=
z" 12.5,fluid is better mixed. At z* = 165.0,more uniformity in the temperature
distribution is achieved and a major portion of the fluid (except near the screw
surface) is at temperatures higher than that at the barrel. This is due to greater
viscous dissipation. At z * = 250.0, i.e. at the exit of the channel, it is clearly seen
that the entire fluid in the x-y plane is at a temperature higher than that at the
barrel. The fluid is now at almost uniform temperature, i.e. thoroughly mixed. It
is interesting to observe that at the exit the isotherms take the shape of stream-
lines (as is evident from vector plots of Fig. 8.27) which was also expected by
Griffith (1962).
 Special Topics II 231
230
Computer Simulationof Flow and Heat Transfer
Table 8.5 Screw Configuration and Input Data
measured, four (including the one at the inlet) are in the metering section while no
measurement is reported for the exit temperature. Out of the six points at which Parameter Value
melt pressures are measured, four are in the metering section including the ones at Diameter, D 63.5 mm (2.5 inch)
the inlet and the exit. In Palit (1974), the experimental points of the pressure and Helix angle, lp 17.66°
temperature profiles are joined by straight lines to show their general behaviour. Metering section channel depth 3.3 mm (0.13 inch)
However, in the comparisons as shown in Fig 8.31 and Fig. 8.32, only relevant Width to height ratio, WIH 16.2
Axial length of the metering 635 mm (25 inch)
experimental data points are indicated.
section, L
Speed, N 100rpm
Power law index, n 0.43
Referencetemperature, To 145°C
Metering section inlet 149 °C
temperature, Tj
Barrel temperature, Tb 150°C
Referencemelt viscosity,Jl.o 6648 N-s/m2
Temperaturecoefficientof 0.009 °el
viscosity, b
Reference strain rate, Yo 1.0 S-I
Flow rate, Q 83 kg/h
Specific heat of the melt, cp 2600 J/kg K
Density of the melt, p 770 kg/m3
Thermalconductivityof the 0.18 W/m K
melt (LDPE), kf
Thermal conductivityof the screw 0.3 W/m K
material (Polyolefin),ks

Table 8.6 Non-dimensional Input Parameters

Parameter Value
Throughput, qv 0.535
Peelet number, Pe 11627.9
Griffith number, G 265.14

For the comparison of the experimental results with numerical data, the mate-
Assumptions Regarding Inlet Melt Temperature In the present analysis,
ial taken is low density polyethylene (LDPE) for which tlle processing tempera-
the inlet melt temperature (assumed uniform in our theoretical model) is taken as
ure is about 150 °e (Palit, 1974). The screw configuration, physical property
149 °e, that is, almost equal to barrel temperature (150 0q. Tj is not taken as
ralues of the melt and the operating conditions of the extruder are given in Table
1.5 for the maximum valve opening meaning that the valve at the exit of the exactly equal to Tb in order to avoid the value of the dimensionless temperature, 8
turning to infinity (see Nomenclature for 8). T!te aforesaid assumption is justified
~xtruder is fully open. The non-dimensional input parameters calculated on the
by the fact that Palit (1974) also constructed his theoretical model of the same
)asis of the data in Table 8.6 are shown in Table 8.6. Q and qv are obtained from
extruder using the developing flow analysis based on the model of Yates (1968)
:he plots of Q vs N and qv vs N respectively corresponding to maximum valve
Jpening as shown in Palit (1974). Furthermore, Palit (1974) reports experiments that assumes that melt being fed through the feed pocket of the extruder and the
inlet melt temperature at the feed pocket being equal to the barrel temperature. It
where the exit valve is either fully open or 20% open. Both situations give rise to
is of interest to note that Palit's model which is based on Yate' s model (1968) is
falling pressure towards the exit and therefore, no comparison could be made for
a quasi two-dimensional one that assumes screw surface temperature to be equal
cases, where pressure is continuously rising towards the exit.
to the barrel temperature which is the "most realistic boundary condition at the
screw surface for a two-dimensional analysis. In other words, it indirectly takes
into account the effect of cross convection. In fact, our three-dimensional model
which does not impose Dirichlet condition at the screw surface also predicts the
screw surface temperature to be almost equal to the barrel temperature which is
attributed to the cross convection effects. Therefore Yate' s model should be
almost as good as the present quasi three-dimensional model.
Calculation of Pressure Profile The pressure profile (Fig. 8.32) has been 2
computed using the experimental metering section inlet pressure (22.6 MN/m )
shown in Palit (1974) as the starting value and the pressure gradient dp/dz at the
inlet calculated earlier numerically.
Comparison of the Present Numerical Predictions with Experimental
Data Figure 8.31 and Fig. 8.32 respectively shows the numerical prediction of
the screw surface temperature and the pressure distribution along the axial length
of the screw expressed in dimensionless forms as against corresponding experi-
mental results of Palit (1974). The predicted dimensional pressure at the exit of2
the extruder is 12.14 MN/m2 as against experimental value of 12.06 MN/m
(Fig. 8.32). Comparison with other two intermediate experimental data points
also show a very good match. The slight deviations may be primarily due to the
assumption regarding the inlet melt temperature and possibly, due to the fact that
234 Computer Simulation of Flow and Heat Transfer

screw body has been considered. An efficient finite difference technique has been
used to obtain the flow, thermal and pressure fields in the screw channel. For the
development of the model, the screw configuration and input data of Karwe and
Jaluria (1990) have been used. The main achievement of this work is in the inclu-
sion of the cross-convection effects (that is, convections normal to the screw
flights as well as to the channel base) in modelling the extrusion process thus
leading to realistic predictions of screw surface temperature and pressure distri-
bution in the downchannel direction. In the absence of experimental data for
Karwe and Jaluria (1990), the numerical results for another extruder of2.5 inch
(63.5 mm) diameter processing LDPE have also been validated satisfactorily with
the corresponding experimental data of Palit (1974). Furthermore, an excellent
comparison of the dimensionless downstream pressure gradient at the end of the
metering section predicted by a quasi two-dimensional computer model of Palit
(1974) with that predicted by the present quasi three-dimensional model also con-
firms the accuracy of the present work.
Nomenclature
b temperature coefficient of viscosity
c specific heat of the fluid
D barrel inner diameter
G Griffith number = 71 V2b/[kj (Tb - Ti)]
H height of the screw channel
Hd depth (below the screw surface) at which the coordinate system is fixed
kj thermal conductivity ofthe fluid
ks thermal conductivity of the screw material
Z axial distance in the metering section ( = z sin cf»
Z* l/L
L axial length of metering section
N screw speed, rpm
n power law index
P pressure
p reference pressure [= 71 (Vb/H)]
Pay average pressure in the x-y plane at each z-location
Pe Peclet number (= Vbz H/a)
Q total volumetric flow rate
qv dimensionless volumetric flow rate or throughput parameter
T temperature
Tb barrel temperature
Ti temperature at the inlet of the heating zone
To reference temperature
u velocity component in x-direction
Vb tangential velocity ofthe barrel (= Tr DN/60)
Vbx component of Vb along x [= Vb sin cf>]
Vbz component of Vb along z [ = Vb cos cf>]
 Special Topics II 237
..:0:>0 computer Simulation of Flow and Heat Transfer

is forced through a die at the extruder discharge where it expands rapidly with reaction. Morgan et al. (1978) have developed a theoretical model which de-
some loss in moisture because of a rapid decrease in pressure. After expanison scribes the effect of temperature-time history, temperature, strain rate, and mois-
and cooling/drying, the extrudate develops a rigid structure and maintains a po- ture content on the apparent viscosity of defatted soy flour dough.
rous texture. It is important to note that flashing or boiling of moisture does not Continued work is required to develop and test the most suitable form of an
occur within the extruder because the pressure exceeds the vapour pressure of apparent viscosity model that adequately, yet simply, describes the important
water at the extrusion temperature (Harper, 1978). changes in viscosity that occur with time during cooking of food dough. Once a
model exists, it will be theoretically possible to simulate the entire cooking extru-
Typical Extrusion Cooked Products Examples of extrusion cooked prod- sion process leading to an increased understanding of the principal variables and
ucts are precooked and modified starches, RTE (ready-to-eat) cereals like corn their influence on extrusion rate, energy input and finished product characteris-
flakes, snack foods, breading substitutes, beverage bases, soft-moist and dry pet tics (Harper, 1978).
foods, processed meat, full-fat and defatted soy flour, macaroni, soup and gravy
bases, and confections like chocolates (Harper, 1978). 8.5 HEAT TRANSFER IN METAL AND ALLOY
The Fundamental Differences Between Polymer and Food SOLIDIFICATION
Processing Although the extrusions of food and thermoplastics closely resem-
ble each other, they have some very sharp differences. One ofthese is the revers- 8.5.1 Introduction
ible melting and solidification of a thermoplastic, as compared to irreversible Metal and alloy solidification find important applications in industry. MetaValloy
cooking reactions for food materials. In cooking, the two most important reac- casting is an industry which involves high annual investments all over the world.
tions, namely protein denaturation and starch gelatinization results in increase of An improvement in the quality of the castings and reducing the number of rejects
dough viscosity as the reactions proceed. In contrast, the apparent viscosity of a would result in a large amount of saving for the casting industry.
thermoplastic will usually decrease as the temperature increases (Fong, 1978). In contrast with the earlier cut -and-try methods, today attention is focussed on
Remsen and Clark (1978) found that for a 25% soy flour suspension, the ap- the scientific aspects of design and production of castings such as the studies of
parent viscosity of the system did decrease initially upon heating until around the heat transfer mechanisms by which castings solidify and the rates at which
160 of (71 0c), then the viscosity started to increase because of the denaturation solidification takes place. Modern researches in this direction are concerned with
reaction (here the starch gelatinization reaction is unimportant). Thus during the the design and placement of risers for castings using computer simulation.
extrusion process, the dough first will encounter "melting" with a decrease in The present discussion deals with (i) one-dimensional pure metal solidifica-
viscosity, then at a certain point down the channel where the dough attains a tion and (ii) two-dimensional alloy solidification. In each case, the method of
certain critical temperature, its apparent viscosity will start to increase demon- solution is briefly mentioned.
strating the "cooking" phenomena. In short, the viscosity of a food system is not
only a function of the strain rate and temperature, but also depend on the compo- 8.5.2 MetalSolidification
sition and the temperature-time history of the process.
Problem Statement (Carnahan et al., 1969) Let us consider a deep region of a
One can consider a plasticating extruder to incorporate all the functions of
solid conveying, melting, mixing and pumping. It is in respect of "melting" and molten metal which is at a uniform initial temperature To. At time t = 0, the tem-
perature of the surface is lowered to Ts, at which it is held constant. The liquid
"cooking" that the analogy between a food and a plasticating extruder stops.
solidifies at a temperature Tj> where Ts < Tf $ To, so that a frozen region grows
Modelling of Food Processing in Single-Screw Extruders From the fore- into the liquid, as shown in Fig. 8.33.
going discussion it should now be clear that to extend the quasi 3-D model of the The objective is to predict the location of the interface and the temperatures in
plastics processing to food extrusion, one has to have a proper viscosity model for regions A and B as functions of time.
the particular food and therein, lies the challenge. Method # 1 Neglecting possible natural convection currents and the small up-
No single model exists describing apparent viscosity as a function of composi- wards motion due to contraction -upon freezing, heat transfer is by conduction and
tion, variability of ingredients, time, temperature and other extrusion parameters. the governing equations are
Earlier researchers have used power-law model to describe isothermal flow
through the metering section for extrusion cooking. Remsen and Clark (1978)
have modelled the time-temperature effects on the apparent viscosity of constarit
moisture soy protein dough which includes activation energy for cooking
8.5.3 Alloy Solidification
Introductory Remarks The physical processes involved in alloy solidifica-
tion are essentially thermal in origin. The present study is aimed at developing an
appropriate thermal model of the solidification of sand-mould casting of an alloy.
It is also intended to calculate the rate of solidification in the casting as a function
of time and predict the occurrence of isolated liquid regions whose shrinkage due
to solidification cannot be compensated by liquid fed from riser. By appropriately
modifying the riser design, it is desired to eliminate casting defects such as shrink-
age cavities with the help ofthe computer model.
Objectives of the Present Study The present study deals with alloy solidifi-
cation within a sand mould with appropriate heat transfer boundary conditions.
The main objective is to determine the riser size and location for producing sound
castings thereby eliminating the costly trial and error method.
The FEM solution technique is used to solve the 2-D transient heat transfer
problem. From the predictions of the analysis, the solidification process within
both the casting and the riser can be visualized by plotting the isotherms at differ-
ent times. From these patterns .the isolated liquid regions which cannot be fed for
shrinkage compensation can be observed. If isolated regions are identified then
the size of the shrinkage cavities can also be predicted. Once the occurrence of
isolated regions is identified, it can be eliminated by changing the size or location
of the riser(s) thereby achieving sound castings. Thus, an optimal riser design
through detailed heat transfer analysis for an alloy casting can be obtained.
Mechanism of Alloy Solidification Before the mathematical modelling of
the problem of alloy solidification is discussed, it is essential to know the salient
features of the process of alloy solidification.
240 Computer Simulation of Flow and Heat Transfer Special Topics II 241

An alloy, unlike a pure metal, does not have a sharply defined freezing tem- as convective heat transfer, specified temperature and in some cases radiative
perature. Instead, it takes place over a range of temperatures. During the process, heat transfer. As a result of these boundary conditions and heat dissipation from
the solid phases separating out at different temperatures possess varying compo- the melt, temperature gradients arise within the melt and the sand mould which
sitions. Due to all these facts, the direction of crystal growth in an alloy depends leads to solidification of the liquid. As the solidification progresses, three distinct
on various factors, such as regions are observed, namely (i) liquid, (ii) mushy, and (iii) solid.
(i) the composition gradient within the casting, These zones have different thermophysical properties.
(ii) the variation of solidus temperature with composition, and
(iii) the thermal gradients within the mould. Boundary conditions at regions 1-10
Consider the phase diagram of a solid solution shown in Fig. 8.35. Let the 1, 9, lO-Constant temperature (T amb)
liquid alloy have the composition Co (of B inA). Also let Tfbe the freezing point 2, 3, 4, 6, 7, 8 - Radiation/convection
of pure metalA, and T/ and Ts' respectively the liquidus and solidus temperatures
(the upper and lower limits for the solidification range) for the alloy composition
Co. As the liquid alloy is cooled down to temperature T/, solid starts to separate
out. The concentration of B in these solids is only C1 as is evident from Fig. 8.35.
As a result, the concentration of B in the liquid, near the solid-liquid interface,
increases to a value more than Co. The solidus temperatures corresponding to the
compositions atP and Q are Tp and TQ' respectively.

Fig. 8.35 Equilibrium phase diagram of A-B alloy

Thermal Model for Alloy Solidification A transient 2-D heat flow model is
developed for alloy solidification with a mushy zone in which the latent heat is
absorbed over a wide range of temperatures. The equations and solution method-
ology developed are general enough to be applicable to a variety of other solidifi-
cation processing problems. A 2-D situation has been modelled for the sake of
simplicity .
The geometry of the sand mould, casting and the sand-base on which the heat
transfer analysis has been performed is shown in Fig. 8.36. The boundaries of the
solution domain are subjected to different modes of heat transfer conditions such
 Special Topics II 243

For castings of various shapes placed inside sand moulds, the heat transfer
analysis is carried out numerically within the entire region of the casting, mould
and sand-base shown in Fig. 8.36. After obtaining the temperature field as a func-
tion of time from such an analysis it is linked to the optimal riser design.
Method of Solution The solution procedure that has been adopted for the
problem under study is the Finite Element Method (FEM) based on Galerkin's
weighted residual approach (See Appendix 'A' for details regarding this method).
Some of the most important advantages are: (a) the easy handling of complex
geometrical shape of the computational domain; (b) the general manner in which
the boundary conditions are introduced; (c) successful representation of compli-
cated material property variations that are difficult to incorporate in other nu-
merical methods.
In view of the afore-mentioned advantages, FEM has been employed to solve
the present problem~
Figure 8.37 shows the finite element mesh for the solution domain which is
actually the half of the physical domain (shown in Fig. 8.36) because of symme-
try. In the present study, 8-noded iso-parametric elements have been used. The
term "iso-parametric" implies that both the element geometry and the variation in
the variable across the element are described utili sing the same type of polyno-
mial function. Finite difference technique has been adopted to approximate the
time derivatives.
Since the pure implicit scheme is the most stable one and has no restriction so
far as the time step is concerned, it has been used in the present study.
Sample Results and Discussion For the sake of brevity, only a few sample
results for aluminium alloy and 0.15% carbon steel bar castings have been shown
and discussed here and the details of optimal riser design have been dropped.
Cylindrical risers are used in all the cases as it is common practice in metal cast-
ing industries. The results have been taken from the M. Tech. thesis of Roy
Chowdhury (1988).
In Fig. 8.38 the process of solidification of a bar casting of aluminium alloy
has been shown. The casting is 85 mm long and 20 mm thick and the riser is
30 mm in dia and 40 mm in height. It is observed that the solidification begins
from the ends of the bar and slowly proceeds towards the riser. This shows that
the thermal gradient is appropriate and leads to directional solidification towards
the source of feed metal which is of utmost importance from the point of view of
the soundness of the casting.
In Fig. 8.38 (a) the isotherm of 625°C is plotted. It shows how after a small
time mushy zone is formed within the casting. The riser retains liquid at this time
which shows that the riser can feed liquid metal to the casting to make up for the
shrinkage which occurs as the casting solidifies. Figure 8.38(b) shows the end of
solidification. The casting has/solidified completely. The isotherm of 525°C is
shown. The mushy zone is within the riser which is desirable. It is clear from the
figure that the riser dimension chosen in this case would give a sound casting
244 Computer Simulation of Flow and Heat Transfer
though this is not optimum riser size. Riser placement is appropriate as the length
of the bar is within the maximum feeding distance.

Fig.8.37 Finite element mesh for the 2-D solution domain

Figures 8.39 and 8.40 show the difference in the solidification rate for bar
castings of steel. In the former case only convective heat transfer loss from bound-
ary is considered while in the latter heat loss due to both convection and radiation
from the exposed surface of the mould is taken into account. It is observed that the
solidification rate is higher due to the enhanced heat transfer when radiation along
with convection is taken into account.
Conclusion
(i) The present problem can be easily extended to more practical situations of
three dimensional geometry.
(ii) The phenomenon of natural convection within the liquid metal can be
looked into and the enhancement of heat transfer due to this phenomenon
can be precisely modelled.
 Special Topics II 247

8.6.2 Need for Thermal Modelling

Detecting thermal problems after component assembly is costly. As a result, PCB
(printed circuit board) thermal analysis which predicts the operating tempera-
tures of the components on the circuit board is found ~o be greatly useful. With the
aid of the numerical model and the results based on it, the designer can redistrib-
ute the components to eliminate localized hot spots, which can lower the average
junction temperature and thus improve reliability (Weiss et aI., 1989).

8.6.3 The Model Problem

The model problem considered here is derived from the literature on cooling of
electronic equipment (Mishra et aI., 1995). The physical situation chosen for
analysis in the present study is depicted
in Fig. 8.41. It is a two-dimensional
idealization of a typical electronic cir-
cuit board array (Kelkar and
. Chowdhury, 1993). The cooling of the
Fig.8.40 Solidification of a bar casting of 0.15% carbon steel for electronic components is taking place
convective and radiative heat loss from the boundary here due to both forced and natural
convection, i.e. mixed convection.
8.6 COOLING OF ELECTRONICEQUIPMENTS Since the channels consist of repetitive
geometric modul.es, the characteristics
8.6.1 Introduction of the flow in successive modules be-
come identical and therefore, analysis
Thermal control of electronic equipments is a challenging task that a heat transfer of an isolated heater will be sufficient.
engineer faces today. This is because with dramatic rise in the density of compo- Figure 8.42 shows such a situation. The
nents on printed circuit boards the power dissipation requirements have jumped flow is assumed to be full developed.
drastically from 5 to about 5 x 104 W1m2 (Weiss et al., 1989). The main objective In the study under discussion, the
of thermal control of electronic equipments is to maintain constant components heated block is of length H and protru-
temperature equal to or below the manufacturer's maximum specified service sion B and is placed in a vertical channel of opening H (Fig. 8.42). The surface of
temperature, typically between 85°C to 100 °C (Peterson and Ortega, 1990). the heater is taken to be at a constant temperature, and the incoming flow in the
When thermal effects go undetected and unchecked, the failure rate of electronic channel is taken to be cold. The walls of the channel are assumed to be thin and
components and assemblies doubles for every increase of 10 to 20°C injunction insulated. Calculations have been carried out at a fixed Prandtl number of 0.7,
temperature (Weiss et at, 1989). with air as the coolant. The method of solution is "operator-splitting algorithm"
Many techniques have been used to cool electronic components. They include described in Chapter 7.
traditional means of natural and mixed gaseous and liquid convection as well as In the following articles, the problem formulation and the method of solution
conduction and radiation or combinations thereof. Some of the other methods are have been presented but results have been given only in synopsis form for the
direct cooling, without phase change, such as that in the Cray-2 supercomputer or sake of brevity. The interested readers are encouraged to refer to Mishra et aI.,
with phase change such as the liquid-encapsulated module where the electronic (1995) for detailed results and their interpretation.
package is submerged in a liquid pool, and indirect cooling schemes, such as the
IBM thermal-conduction module (TCM) in the IBM 3081/3090 processor or the Problem Formulation Equations governing flow and heat transfer are ex-
liquid-cooled module (LCM) in the NECSX. Other more advanced methods pressed in terms of stream function 1/1, vorticity ~, and temperature T in dimen-
include miniature thermosyphons or free~falling liquid films (Peterson and sionless form as follows (Roache, 1972):
Ortega, 1990).
 252 .Computer Simulation of Flow and Heat Transfer Special Topics II 253

SUMMARY 6. Fenner, R T, "Developments in the Analysis of Steady Screw Extrusion of

Polymers", Polymer, Vol. 18, 1977, pp 617--635.
In Chapters 4-7 the application of computational techniques to basic heat transfer 7. Fong, DSC, Experimental Study of Extrusion-Cooking of Defatted Soy Flour,
and fluid flow processes was discussed. Except in Chapter 4 (Sec. 4.27.2) where M.S. Thesis, Chemical Engineering, Virginia Polytechnic I.nstitute and State
.handling of non-linearity in the boundary condition was demonstrated through a University, BaIcksburg, VA, USA, 1978.
plane wall conduction problem having one boundary exposed to radiation heat 8. Ghoshdastidar, P S, C A Rhodes, and D I Orloff, "Heat Transfer in a Rotary Kiln
during Incineration of Solid Waste", 23rd ASME National Heat Transfer
flux, no treatment of numerical methods for radiation heat transfer was given. In
Conference, Denver, Colorado, August 4-7, AS ME paper No. 85-H5-86, 1985.
this chapter, numerical methods for calculation of radiation in enclosure are dealt
9. Ghoshdastidar, P Sand YST Raju, ''Transient Combined Mixed Convection and
with for rotary kiln problems. The absorption factor method is described in detail Radiation from a Vertical Aluminium Fin", Journal of Energy, Heat and Mass
in AppendixD. The combined modes of heat transfer such as convection and Transfer, Vol. 14, 1992, pp 45-53.
radiation are demonstrated in detail for the vertical fin problem and the alloy 10. Ghoshdastidar, P S and V K Anandan Unni, "Heat Transfer in the Non-Reacting
solidification problem and mentioned briefly for the case of outer wall of the Zone of a Cement Rotary Kiln", ASME Journal of Engineering for Industry, Vol.
rotary cement kiln. Another example of combined mode is shown in the problem 118 (February issue), 1996, pp 169-172.
of the cooling of electronic equipments where both free and forced convection 11. Gilbert, W, "Heat Transmission in Rotary Kilns", Cement and Cement
(mixed convection) are occurring. Conjugate heat transfer is demonstrated in the Manufacture, Vol. 9,1936, pp 115-128, pp. 139-154 pp 207-220.
verticalfin as well as in the screw extruder problems. The numerical treatment of 12. Griffith, R M, "Fully Developed Flow in Screw Extruders", Ind. Eng. Chem.
non-Newtonian.flow and heat transfer is shown in the case of polymer and food Fundam., Vol. 1, No.3, 1962, pp 181-187.
13. Harper, J M, Extrusion of Foods, CRC Press, 1978:
processing in single-screw extruders. The computational methods for handling
14. Helmrich Hand K Schtigerl, "Rotary Kiln Reactors in Chemical Industry", Ger.
phase change in metal and alloy soldifications are described. The various prob-
Chem. Engg., Vol. 3, 1980, pp 194-202.
lems discussed here have been solved by different numerical methods such as 15. Holman, J P, Heat Transfer, 5th International Student Edition, McGraw-Hill,
finite-difference, finite-element and control volume methods. The problems ana- 1981.
lysed in this chapter range from enhanced heat transfer, manufacturing processes 16. Hottel, HC, Radiant Heat Transmission, Chap. 4,. in W.H. McAdams (Ed.),
for metal, alloy, polymer and food, applications to drying and solid waste treat- "Heat Transmission", 3rd Ed., McGraw-Hill Book Company, New York, 1954.
ment in rotary kilns to cooling of electronic components. A detailed discussion of 17. Jacob,.M, Heat Transfer, Vol. II, John Wiley & Sons, Inc., New York, 1957, pp
results follows each problem formulation to emphasize the need for design and 19-21.
optimization of complex industrial systems involving thermal processes. 18. Karwe, M V andY Jaluria "Numerical Simulation of Fluid Flow and Heat Transfer
in a Single-Screw Extruder for non-Newtonian Fluids", Numerical Heat Transfer,
Part A, Vol. 17, 1990, pp 167-190.
REFERENCES 19. Marshall, D I, I Klein, and R H Uhi, SPE J., Vol. 21, 1985, p. 1192.
20. Mishra, Debasish, K MuraIidhar, and P S Ghoshdastidar, "Computation of Flow
1. Carnahan, Brice, H A Luther, and J 0 Wilkes, Applied Numerical Methods, John and Heat Transfer Around a Vertical Discrete Protruding Heater Using an
Wiley & Sons, New York, 1969. Operator-Splitting Algorithm':, Numerical Heat Transfer, Part A, Vol. 28, 1995,
2. Carslaw, H Sand J C Jaeger, Conduction of Heat in Solids, 2nd Ed., Oxford pp 103-119.
University Press, New York, 1959. 21. Morgan, R G, D A Sutter and V E Sweat, "Modelling the effects of temperature-
3. Chowdhury, Dipankar, and Kanchan M Kelkar, "Numerical Prediction of time history, temperature, shear rate and moisture on apparent viscosity of defatted
Periodically· Fully Developed Natural Convection in a Vertical Channel with soy flour dough", Paper No. 79-6002, Armerican Society of Agricultural
Surface Mounted Heat Generating Blocks", Int. J. Heat Mass Transfer, Vol. 36, Engineers, St. Joseph, MI, 1978.
No.5, 1993, pp 1133-1145. 22. Palit, K, Some Experiments on Polymer Melt Flow in a Plasticating Extruder,
4. Das, Manab K, A Detailed Numerical Study of Thermal Transport Processes in Polymer Processing Research Report No.1, Imperial College, Mechanical
the Metering Section of a Single-Screw Plasticating Extruder, Ph D Thesis, Engineering Department, U K, 1974.
Mechanical Engineering, lIT Kanpur, India, 1993. 23. Patankar, S V, Numerical Heat Transfer and Fluid Flow, Hemisphere,
5. Das, Manab K and P S Ghoshdastidar, "A Quasi Three-Dimensional Computer Washington, DC, 1980.
Model of Flow and Conjugate Heat Transfer in the Metering Section of a Single- 24. Patankar, S V and D B Spalding, "A Calculation Procedure for Heat, Mass and
Screw Plasticating Extruder", Proc. 1994 ASME/AIAA Thermophysics and Heat Momentum Transfer in Three-dimensional Parabolic Flows", Int. J. Heat and
Transfer Conference, Colorado Springs, Colorado, June 20-23, ASME HTD- V01. Mass Transfer, Vol. 15, 1972, pp 1787-1806.
275, 1994, pp 123-133.
254 Computer Simulation of Flow and Heat Transfer

25. Peterson, G P and A Ortega, "Thermal Control of Electronic Equipment and

Devices", Advances in Heat Transfer (Ed. Hartnett, J P and Irvine, Thomas F, Jr),
Vol. 20, 1990, pp 181-314, Academic Press Inc.
26. Raithby, G D and G E Schneider, "Numerical Solution of Problems in
Incompressible Fluid Flow: Treatment of the Velocity-Pressure Coupling",
Numerical Heat Transfer, Vol. 2, 1979, pp 417-440.
27. Remsen, CHand J P Clark, "A Viscosity Model for a Cooking Dough", J. Food
Process Engineering, Vbl. 2, No.1, 1978, pp 39-64.
28. Roache, P J, Computational Fluid Dynamics, Chap. 5, Hermosa Publishers,
Albuquerque, New Mexico, 1972.
29. Roy Chowdhury, Saibal, A Finite Element Based Thermal Analysis of Optimal
Riser Design for Alloy Castings, M. Tech. Thesis, MechanicalEngineering, lIT
Kanpur, India, 1988.
30. Sass, A, "Simulation of the Heat Transfer Phenomena in a Rotary Kiln", Ind. Eng.
Chem. Process Design and Development, Vol. 6, No.4, 1967; pp 532-535.
31. Schlichting, H, Boundary Layer Theory, McGraw-Hill, New York, 1979.
32. Spang, H A, "A Dynamic Model of a Cement Kiln", Automatica, Vol. 8, 1972,
pp 309-323.
33. Sunden, B, "Conjugate Mixed Convection Heat Transfer from a Vertical A.I THE ESSENCE OF THE FINITE ELEMENT
Rectangular Fin", Int. Comm. Heat and Mass Transfer, Vol. 10, 1983, pp MEmODS
267-276.
34. Suryanarayana, N V, J F Lyon, and N K Kim, "Heat Shield for High Temperature Finite element methods (FEM) basically seek solutions at discrete spatial regiol1.s
Kiln", Ind. Eng. Chern. Process Design and Development, Vol. 25, 1986,
pp 843-849. (called elements) that fill the solution domain by assuming that the govemil1.g
differential equations apply to the continuum within each element. Many convel1._
35. Van Doormaal, J P and G D Raithby, "Enhancements of the SIMPLE Method for
Predicting Incompressible Fluid Flows", Numerical Heat Transfer, Vol. 7, 1984, ient shapes for the element are available, such as triangles and quadrilaterals in
pp. 147-163. two dimensions and tetrahedra, pentahedra, and hexahedra in three dimensions.
36. Weber, P, Heat Transfer in Rotary Kilns, Bauverlag Gmbh, Wiesbaden-Berlin, FEM is based on integral minimization principle and provides piecewise (Or
1963. regional) approximations to the governing equations.
37. Weiss, Jonathan, Patrik Fortner, Bruce Pearson, Keith Watson, Tom Monroc, The finite element method was originally conceived by engineers in the 1950's
"Modelling Air Flow in Electronic Packages", ASME Mechanical Engineering, to analyse aircraft structural systems using the emerging digital computer. As tl)e
September, 1989, pp 56-58. FEM in solid mechanics gained ground, the concept of a "force balance" Wils
38. Yates, B, Temperature Development in Single Screw Extruders, Ph D Thesis, - replaced by theory founded in the variational calculus and classical Rayleigh_
University of Cambridge, UK, 1968.
Ritz methods (Baker and Pepper, 1991). The introduction of the variationill
calculus based FEM and its ready acceptance in solving linear heat transfer and
fluid mechanics problems such as conduction and creeping flow was dUe to
relative ease with which a problem with complicated boundary shapes cab be
modelled, especially when compared with finite-difference methods. Anoth~r
advantage is the handling of thermal stress problems for which temperature calC\!_
lations are required. Therefore, using the same "elements" the heat transfer and
the stress engineers can solve different aspects of the same problem (Myers,
1971). However, this method cannot be extended to convective heat transfer and
high Reynolds number fluid flow problems which are essentially nonlinear in
nature. This is because the fluid convection term couples velocity to field
gradients (the field can be velocity in case of momentum trlYlsportor temperatUte
in case of energy transport) which does not admit creation of a variation~l
principle as required by the Rayleigh-Ritz theory.
A.9 ACCURACY OF SOLUTION

Comparison of the FEM solution with the exact solution (Eq. (A.51)) at the nodal
points 1,2, 3 reveals exact matching between the two. However, this kind of
exact agreement occurs only when a problem is very simple as the present one. In
other words, it can be said that the use of piecewise linear trial functions has
correctly interpolated a quadratic function, indicating the good performance ex-
pected in the general case. It must be noted, however, that any numerical solution
is an approximation, and hence will not agree everywhere with the exact solution.
The accuracy of the approximate solution depends on the location of evaluation.

A.tO SUMMARY

The finite element method based on Galerkin's weighted residual approach is

described in detail for a one-dimensional steady state conduction problem. For
more details regarding multi-dimensional steady state problems and unsteady
transport with fluid motion, readers are referred to Baker and Pepper (1991).
280 Computer Simulation of Flow and Heat Transfer

A.II CONTROL VOLUMEFORMULATION

(ALSO KNOWN AS FINITE VOLUME
MEmOD)

A.II.I Basic Concept

In the control volume fonnulation, the computational domain is divided into a
number of nonoverlapping control volumes such that there is one control volume
surrounding each grid point. The governing differential equation (GDE) is inte-
grated over each control volume. Piecewise profiles expressing the variation of
the unknown function (such as temperature, 1) between the grid points are used to
evaluate the required integrals. This results in the discretization equation consist-
ing of the value of the unknown function for a cluster of grid or nodal points.

A.II.2 One-Dimensional Steady-State Heat

Conduction Problem
Let us consider steady state I-D heat conduction in a medium of conductivity 'k'
(uniform or non-unifonn) in which heat is generated at the rate of q'" W/m3•
The governing differential equation is:
A.ll.6 Summary
The control volume formulation for the 1-0 steady state conduction problems is
described in detail. The extension to 2-D and 3-D situations are discussed briefly.
For transient conduction, flow and convection-diffusion problems, readers are
referred to Patankar (1980).

A.12 CLOSURE

In this chapter three discretization methods other than the finite-difference

method, such as variational calculus based FEM, Galerkin's weighted residual
approach based FEM and control volume method are discussed in detail with
respect to one-dimensional steady state conduction problems. Interested readers
are encouraged to advance their knowledge in these areas by going through the
relevant books given in the references.

REFERENCES

1. Baker, A J and D W Pepper,Finite Elements 1-2-3, McGraw-Hill Inc., New York,

1991.
2. Galerkin, B G, "Series Occurring in Some Problems of Elastic Stability of Rods
and Plates", Engrg. Bull., Vol. 19, 1915, pp. 897-908.
3. Myers, Glen E, Analytical Methods in Conduction Heat Transfer, McGraw-Hill,
New York, 1971.
4. Patankar, S V, Numerical Heat Transfer and Fluid Flow, Hemisphere,
Washington, DC, 1980.
c.t SUBROUTINE TDMA
c SUBROUTINE FOR SOLVING A SYSTEM OF LINEAR
c SIMULTANEOUS ALGEBRAIC EQUATIONS HAVING
c A TRIDIAGONAL COEFFICIENT MATRIX.
c THE EQUATIONS ARE NUMBERED FROM I THROUGH N,
c AND THEIR SUB-DIAGONAL, DIAGONAL AND SUPER-
c DIAGONAL
c ELEMENTS ARE STORED IN THE ARRAYS A, B, AND C. THE
c RIGHT-HAND COLUMN VECTOR ELEMENTS ARE STORED
c IN THE ARRAY D. THE COMPUTED SOLUTION VECTOR XCI) ".
c X(N)
c IS STORED IN THE ARRAY X.
SUBROUTINE TDMA (I, N, A, B, C, D, X)
DIMENSION A(1), B(1), C(1), D(1), X(1), BETA (101),
GAMMA (101)
c *** COMPUTE INTERMEDIATE ARRAYS BET A AND GAMMA **
BETA (I) = B(I)
GAMMA (I) = D(I)/BET A(I)
11=1+1
DO 1 J = 11, N
BETA (J) = B(J)-A (J)*C(J-1)/BETA(J-1)
1 GAMMA(J) = (D (J) - A(J) * GAMMA (J-1)/BETA (J)
c ***** COMPUTE THE SOLUTION VECTOR X *
X(N) = GAMMA (N)
N1 = N-I
DO 2 K = 1, N 1
J = N-K
 292 Computer Simulation of Flow and Heat Transfer

2 X (J) = GAMMA (J) - C(J) * X (J + 1) /BETA (J)

RETURN
END

C.2 DEMONSTRATION PROGRAM SHOWING

APPliCATION OF SUBROUTINE liMA
The application of subroutine TDMA is demonstrated through the solution of
Eq. (4.40). In the demonstration program A, B, C, D are the coefficient vectors
defined in Eq. (4.23). Tis the vector of the unknowns. In the subroutine TDMA,
BET A and GAMMA are vectors of intermediate coefficients Pi and 'ii, respec-
tively defined in Eq. (4.28). The level of programming in the demonstration pro-
gram has been kept intentionally simple to make the readers appreciate how the
arrays A, B, C, D are defined. The output is given below the main programt
c DEMONSTRA nON PROGRAM SHOWING APPLICA nON OF
c SUBROUTINE TDMA
DIMENSION A(21), B(21), C(21), D(21), T(21)
OPEN (UNIT = 76, FILE = 'OUTP')
A(2) = - 1.
A(3) = - 1.
A(4) =-2.
B(1) = 2.25
B(2) = 2.25
B(3) = 2.25
B(4) = 2.25
C(l) = - 1.
C(2) = - 1.
C(3) = - 1.
D(1) = 1.
D(2) = O.
D(3) = O.
D(4) = O.
CALL TDMA (1, 4, A, B, C, D, T)
WRITE (76, 2) (T (1),1 = 1,4)
2 FORMAT (2X, 4F8. 3)
STOP
END
.629 .415 .305 .271

REFERENCE

1. Carnahan, Brice, H A Luther, and James,a. Wilkes, Applied Numerical Methods,

John Wiley & Sons, New York, 1969.
 The Absorption Factor Method 295

4. Gebhart also proved a reciprocity relationship between the absorption

factors. For diffuse radiation and reflection:
A-I = C·J B··JI A-J
C·I B··
I}

5. The absorption factors depend only on the geometry, through sh~pe factors
and the properties C and P of the surfaces comprising the enclosure. They
do not depend on the surface temperature and the heat input. This is
particularly advantageous because Bjj values are obtained independent of
the surface temperature and the heat input. Therefore, in many complicated
problems of thermal design of industrial systems, the absorption factor
Fig. D.I An enclosure of N surfaces method is more suitable compared to other methods since the numerical
solution for Bij is carried out conveniently.
Knowing the shape factors and the values of the reflectivity and the emissivity,
the values of the Bij's can be determined numerically. The follo~ng equations
can be written for each of the N surfaces. REFERENCE
BIj = Flj £j + Fll PI BIj + Fiz Pz BZj + ... + FIN PN BNj
BZj = FZj Cj + FZI PI BIj + Fzz Pz BZj + ... + Fw PNBNj 1. Gebhart, B., Heat Transfer, 2nd Ed., McGraw-Hill Book Co., New York, 1971.

BNj = FNj Cj + FNl PI Blj + FNZPZ B2j + ... + FNN PNBNj (D.2)
=
wherecj, Pi' i 1... N represent emissivity and reflectivity oftheith surface. F y's
are the shape factors.
These N linear equations can be transposed and rearranged to obtain:
(Fll PI -l)BIj + FIZ Pz BZj + ... + FIN PNBNj + FIj=0 Cj

F21 PI BIj + (Fzz PZ - l)BZj + .. , + Fw PN BNj + FZj E.i= 0

FNl PI Blj + FNZPZBZj + ... (FNN PN-l) BNj + FNj Cj = 0 (D.3) .
Equation (D.3) can be solved for the N unknowns Blj' BZj' BNj' Equation (D.3)
is valid for any choice of j, that is, 1, 2, ... N. For small N, the direct methods for
linear simultaneous equations such as Gaussian Elimination may be employed.
When N is large, iterative method such as Gauss~Seidel should be used.
With the values of Bij known, the radiant heat loss or gain qj at each wall
element is calculated using Eq. (D. 1).
Some special points of importance should be noted:
1. The distinction between the Bij and F ij is as follows. The shape factor
represents the fraction of the radiant energy leaving the area Ai and
arriving directly at the area Aj while the absorption factor represents the
fraction of the radiant energy leaving the areaAj and absorbed by the area
Aj following multiple reflections before it finally arrives at the areaAj.
=
2. In the special case of black surfaces: B ij F ij'
3. Since the energy emitted by each surface is absorbed by the collection of
N surfaces which form the enclosure,
 Index 297

CTCS 71 Flow over flat plate 140·

Cylindrical geometry 90 Flow over a heated flat plate 162
Flow of liquid metals 163
Derivatives 7 Flow over 2-D wedge 140
Diagonal dominance 45, 187 Forward difference 21,23,24,25
Difference equation 1 Fourier, Joseph, J.B. 3
Differential equation 1,7,8 Free convection 170
Diffuse-gray surfaces 293 Free convection from a vertical
Diffusion equation 11 plate 170
Dirichlet condition 13 Free convection in enclosure 174
Discretization error· 29 FTCS 71
Discretization, methods of 7
Absorption factor method 293 Cauchy condition 14 DufortcFrankel method 81 Gaussian elimination 1, 4, 40, 48, 52,
Accuracy 21,25,27,28,40,70,279 Cement rotary kiln 206 59
AD! 82,84,86,89 Central difference 21,22,23,24,25,27 Eigen-values 74, 75 Gauss-Seidel iterative method 40,42,
Advantages of computer simultation 2 CFD 7 Eigen-vectors 74 44,48,51,52,59
Age of earth 3 Characteristics 15 Electronic equipments, cooling of 246 G-E 40
Algebraic equations 1,4 Checkboard 133 Elliptic PDE 11 Gebhart, B. 293
Algorithm: Closed domain 11 ENIAC 7 Grashof number 201
Operator-S plitting(OS) 183, 184, Combined modes of heat transfer 190 Equation: G-S 42,43,45,89
189,249 Compatibility conditions 97,98 Biharmonic 126 Global diffusion matrix 277
Thomas 40,41,42 Composite fin 179 Continuity 120 Global load matrix 277
Tridiagonal Matrix (TDMA) 40,41, Composite media 96 Euler-Lagrange 259 Global matrix 277
42 Computational fluid dynamics 7 Laplace's 11 Governing differential equation,
SIMPLE 137 Computer simulation 1, 8 Momentum 120 nonlinear 98
Alloy solidification 239 Condition at the centre, handling of 90 Navier-Stokes 7, 119, 120, 122, 132 Grid 4,5
Alternating Direction Implicit Conduction heat transfer 32 Poisson's 11, 129, 137 Grid independence test 29,87,88
Method 82 Consistency 80 Pressure-correction 137 Grid points 4
Alternative to upwind scheme 183 Continuity equation 120 Vorticity transport 122
Analytical method 5,6 Control volume formulation 8,280 Wave 12 Heat conduction 11
Application of ADI to composite Convection-conduction parameter, Euler~Lagrange equation 259 Higher accuracy 25
fin 179 CCP 192 Euler(or Explicit) method 66, 70, 71, Hybrid scheme 158
Application of subroutine TDMA 292 Convective (or convection) boundary 72,73,74,75,77,78,79,87,88, Hyperbolic PDE 11
Axisymmetric problems 90 condition 53, 76,265 100, 145, 163, 165, 172
Convection heat transfer 153 Exponentialscheme· 158 Image-point technique 35, 37, 38, 39,
Backward difference 21,24,25 Convection-diffusion 154, 183, 185, 54,93
Back substitution 49 287 False diffusion 157 Implicit 163, 174
Banded coefficient matrix 59 Convection-diffusion (steady, one- False transient 89 Implicit Crank-Nicholson method 166
Band width 59 dimensional) 154 FEM 7,259,265 Implicit Keller-box method 167, 170
Biharmonic equation 126 Convection-diffusion (unsteady, one- Fin 33, 179, 190 Incineration of solid waste 214
Block elimination 170 dimensional) 158 Finite-difference 1,7 Incompressible fluid flow 119
Boundary conditions 11,13 Convection-diffusion (unsteady, two- Finite element 7 lildustrial problems 190
Boundary conditions, nonlinear 3,101 dimensional) 159 Finite element basis 271 Initial condition 12,13 .
Boundary layer 15, 138, 140 Cooling of electronic equipments 246 Finite element method, essence of 255 Initial value problem 15
Boundary layer flow 138, 144 Comer points, handling of 57 Finite element method based on Insulation cO!1dition 14
Boundary value problem 15 Crank-Nicholson method 66, 68, 70. variational ca:Iculus 256 Integral minimization 7,261
BTCS 71 71, 72, 73, 74, 75, 78, 79, 86, 87, 166 Finite volume method 280 Interpolation function 272
Calculus of variations 256 Creeping flow 125 Flops I Inviscid flow 127
 ~,lt Index
Index 299
Irregular geometry 3,106 Peclet number 157
Implicit Keller-box 167, 170 Solution of higher order ODE 290
Irrotational flow 122
Line-by-Line 60,61 Phase change 190 SOR 47, 89
Newton-Raphson 143 Poisson's equation 11 Source-term linearisation 282
Jacobi method 43,44
Pure implicit 69 Poisson equation for pressure 129 Spalding D.B. 132
Relaxation 46 Poisson equation for pressure-correc- Spectral method 8
Kronecker delta 275 Runge-Kutta 288 tion 137 Spherical geometry 95
Shooting 142 Polymer and food.extrusion 126 Stability 68,71,77,78, 173
Lagrange interpolation 272 Polymer and food processing 220
Spectral 8 Stagg~red grid 132, 134
Laplace's equation 11 Polynomial 7
Methods of discretization 7 Stream function 122
Laser cutting 16
Mixed convection 153, 190 Polynomial fitting 37 Stream function-vorticity method 121
Laser surgery 16
Mixed convection and radiation 190 Prandtl, L. 138 Subroutine TDMA 292
L'Hospital's rule 92
Mixed elliptic-parabolic Prandtl number 154,201 Subsonic 16
Linearity 10 122
Modelling 190 Pressure-correction equation 137 SUR 47
Line-by-Line method 60,.61
Momentum equations 120 Pressure-correction equation, Boundary Supercomputers 5
LOCA 2
conditions for 137 Supersonic 16
Natural convection 170 Pressure-gradient 120
MAC 132
Navier-Stokes equations 7, 119, 120, Power-law scheme 158 Taylor series 7,21,81,97
Marker-and~Cell 132 122; 132 Primitive-variables approach 131 TDM 39,59
Matrix:
Neumann condition 14 Problem complexity 4 TDMA 40,41,42,48,49,52,159,187
banded coefficient 59 Pure Implicit Method 69
Newton-Raphson method· . 143 Thermal boundary layer 162, 163
.global 277
Non-axisymmetric problem 93 Thermal radiation 207,293
global diffusion 277 Radiation boundary condition 103
Non-Fourier conduction 16 Thomas algorithm 40,41,42
global load . 277 Radiation in enclosure 190,293
NOnlinear boundary conditions 3,101 Three-dimensional problems 62, 267,
tridiagonal 39 .
Nonlinear governing differential Recursion 41,47 286
tridiagon(ll block 170 Relaxation factor 47
equation 3,98 Three-dimensional transient
Mega 1
Nonlinearlity 98, 120 Relaxation method 46 problems 89
Magaflop 1;5
Non-Newtonian flow 190 Reynolds number 7, 154 Transient one-dimensional problems 63
Metal and alloy solidification 237 Robbins condition 14
Non-uniform grid 25 Treatment of source term 281
Metal solidification 237
Numerical errors 21,28 Rotary kiln, cement 206 Trial function 270
Materials processing 3,16
Numerically induced oscillations 71 Rotary kiln, heat transfer in 205 Tridiagonal block matrix 170
Method:
Numerical method 5,6 Rotary kiln incinerator 214 Tridiagonal matrix 39
Absorption factor 293 Rotational flow 122
Numerical simulation 2 Tridiagonal matrix algorithm 40,41,42
Alternating direction Round-off error 28 Truncation error 22,29
Implicit(ADI) 82 Runge- Kutta methods 288
One-dimensional steady state Two-dimensional steady state 54, 266,
Control volume 8,280 33, 259,
268,280,283 Runge-Kutta methods, essence of 288 286
Crank-Nicholson 66,68, 70, 71, 72, Operator-splitting algorithm 183, 184,
73,74,75,78,79,86,87,166 189,249 Scarborough criterion 44,45· Upwind scheme 156, 157, 159, 161
Dufort-Frankel' 81 Similarity solution 140
Open-ended 11,12 Upstream differencing 156
Euier (or Explicit) 66, 70; 71, 72, 73,
Optimum relaxation factor 47 Shooting method 142 Under-relaxation 46
74,75,77,78,79,87,88,100, Shooting technique 140 Uniforrii grid 21
Optimum step size 29
145, 163, 165, 172 SIMPLE' 132, 135, 137
Order 10
Finite-difference 1,7,21
Oscillations, numerically induced SIMPLE algorithm 137 Variable thermal conductivity 98
Finite element 7,256,268 71
Over-relaxation 46 SIMPLER 138 Variational principle based FEM 7.256
Gaussi.an elimination 1,4,40,48, SIMPLER revised 138 Vertical fin 190
52;59 Simultaneous ordinary differential
PDE 10 Vertical plate, Free convection
Gauss-Seidel iterative 40, 42, 44, equations 289 from 170
Parabolic PDE 11
48,51,52,59 Single-screw plasticating extruder 220
Partial differential equations 10 Von-Mises transformation 147
Implicit Crank-Nicholson 166 Small Reynolds number 125 Vorticity 121
Patankar S.V. 132, 135
300 Index

Vorticity transport equation 122 Weighted residual approach,

Galerkin's 7,268
Wave equation 12 Weight function 270
Weak statement 272 Wiggles 156