Volume 3 No.2                                                                                                                     January 2000


Historical Developments in Computational Fluid Dynamics (CFD)

1. The Subject

The physical aspects of transport phenomena in the macroscale are governed by Newton's laws of motion and the fundamental principles of mass, energy and species conservation. Depending on the nature of the problem and the quantities of interest, these fundamental principles can be expressed in terms of algebraic equations, ordinary/partial differential equations or integral representations. Numerical simulation is by and large the technique of replacing the governing transport equations with algebraic equations and obtaining a final numerical description of the phenomenon in space and/or time domain. Irrespective of the nature of the problem, numerical simulation involves the manipulation and solution of numbers, leaving behind an end product which is also a collection of numers. This is in contrast to the symbolic expression of closed form analytical solution.

The final objective of most engineering investigations, whether they are analytical, experimental or numerical, is to obtain a quantitative description of the problem, in terms of numbers. In this regard, numerical simulation techniques provide readily acceptable and often the most descriptive form of solutions to a variety of transport problems. Numerical simulation of practical problems generally involves the repetitive manipulation of thousands, or even millions, of numbers - a task that is feasible only with the aid of a computer. Advances in simulation techniques and their applications to problems of greater complexity and enormity, are intimately related to the advancement in computer technology. This is exemplified by the fact that one of the strongest forces driving the development of new supercomputers is the increasing demand for speed and storage required by the numerical simulation of the flow problems.

2. Historical Perspective

By the turn of the twentieth century, the development of closed form analytical solutions for field problem had reached a highly mature stage and it was being realized that a large class of problems still remained which were not amenable to exact analytical solution methods. This gave birth to a variety of approximate semi-analytical techniques on the one hand and to the development of numerical solution procedures on the other. The semi-analytical techniques which foundwide use in fluid dynamics research were the perturbation methods [1], the similarity approach [2], and the integral method [3], - all for the viscous boundary layer calculations, - and the method of characteristics [4] for inviscid compressible flow simulations. As regards the numerical techniques for solving field problems, finite difference based methods were the first to develop, because of their straight forward implementation. Although the finite difference formulation is relatively simple, the severe limitation faced in the pre-second world war era was that calculations had to be performed manually. Thus, even linear problems involving Laplacian or Biharmonic operators were solved iteratively by relaxation methods [5,6]. Southwell [7] introduced a relaxation scheme which is highly suitable for hand calculations. In this method, the residuals of the governing equations are calculated at all the grid points of the solution domain and the variable values corresponding to the locations of largest residuals are relaxed first. Until the advent of digital computers, Southwell's method was very popular for solving various heat transfer and fluid flow problems. Another relaxation scheme which found extensive use was the Successive Over Relaxation (SOR) scheme proposed by Frankel [8].

For structural problems involving elastic deformations, Ritz [9] developed a method which involves the approximation of a potential functional (virtual work) in terms of trial functions with undetermined coefficients. The unknown coefficients are evaluated by minimizing the potential functional. A severe limitation of the Ritz method is that the trial functions need to satisfy the boundary conditions of the problem. Courant in 1943 [10] made a significant improvement over Ritz method by discretizing the domain into triangular areas and assuming linear trial functions over each of the triangles. By this ingenious extension, all the trial functions were not required to satisfy the boundary conditions. Incorporating these concepts, the full fledged development of the Finite Element Method was first introduced by Clough [11] in 1960. Since then, the method has made rapid strides for the modeling of structural engineering problems and fluid flow and heat transfer modeling in recent years.

A pioneering work on the uniqueness and existence of numerical solutions to partial differential equations was presented by Courant, Friedrichs and Lewy[12] in 1928. The stability requirement (CFL condition) for hyperbolic partial differential equation was first shown in this work. The stability criteria for Parabolic time-marching problems in Computational Fluid Dynamics were developed by Von Neumann. A detailed discussion of these criteria has been presented by O'Brien et al [13]. With so much of a ground work having been accomplished prior to 1950 on the basic numerical methods, iterative schemes and numerical stability, the progress on numerical simulation was accelerated by leaps and bounds after the discovery of the electronic computers in the late fifties.

One of the earliest flow problems to be attacked with the help of digital computer was the viscous flow simulation at intermediate Reynolds numbers (Re < 1000). Based on the stream function-vorticity formulation of viscous flow problems, Fromm [14] and Fromm and Harlow [15] developed an explicit forward time difference method at Los Alamos. Their method was used by Thoman and Szewcsyk [16] for cross flow over cylinders and Rimon and Cheng [17] for uniform flow over a sphere. Steady state solutions for stream function-vorticity equations were obtained by Hamielec et al [18,19] using the SOR technique. An implicit time-marching procedure for viscous flows was developed by Pearson [20]. This method is based on the Alternating Direction Implicit (ADI) method proposed by Peaceman and Rachford [21] and Douglas and Rachford [22].

Significant development was witnessed in the fifties and sixties towards the solution of inviscid compressible flow equations. Starting with shock-capturing technique of Lax [23] which used the conservative form of governing equations, several methods have been developed. Finite difference schemes, such as the Particle-in-Cell (PIC) have been found to be inherently shock smearing [24]. In 1960, a second order accurate finite difference scheme which reduces the shock smearing effect was proposed by Lax and Wendroff [25]. This scheme later led to the development of the McCormack method [26]. For moving shocks, shock fitting procedures have been proposed and applied to multidimensional supersonic flows over various configurations [27, 28, 29]. Even today, some of these schemes are in extensive use.

Although in the early simulation methods for viscous incompressible flow, vorticity and stream function were the calculated variables, and the late sixties, simulations in terms of primitive variables (velocity components and pressure) began. Pioneering work in this direction was performed by Harlow and Welch [30] and Harlow and Amsden [31] at Los Alamos. These authors introduced explicit transient algorithms such as MAC and SMAC. Chorin [32] in 1968 developed the artificial compressibility method for handling viscous incompressible flows. Adopting some of the concepts proposed in these studies, a successful implicit formulation in terms of primitive variables was developed by Patankar and Spalding [33]. Based on this well known SIMPLE algorithm and its later improvements such as SIMPLER [34] and SIMPLEC [35], a horde of multi-dimensional viscous incompressible flows have been simulated. These implicit methods have an inherent advantage over the explicit algorithms that they have no restrictions on the time step from the point of view of numerical stability.

In the late seventies and eighties considerable interest has been evinced on the techniques for handling flows in arbitrary shaped geometries. Methods for transforming complex geometries into simple ones have been proposed and excellent discussion on these methods are provided in the book by Thompson, Warsi and Mastin [36] and the reviews on this subject [37, 38]. In recent years, Baliga and coworkers [39, 40, 41] have introduced a control volume based finite element method which can handle arbitrary geometric. Independently, weighted residual based finite element algorithms have been developed [42, 43, 44]. Peric [45] and Majumdar et al [46] have proposed the application of conservation laws to non-orthogonal control volumes. Research is still being conduction on these methods and several complex flow situations are being simulated using them.

Although a large volume of research publications have appeared in recent years on numerical flow simulation, the potential for further research is expanding at an ever increasing rate. In the decades to come, it appears that many more powerful algorithms will be evolved and several complex flow/heat transfer problems will be successfully simulated. The scope and potential for applying numerical simulation in real life problems is described briefly in the next section.

3. Role of Numerical Simulation in Modern Technological Environment

The role of numerical simulation in engineering is so vital that it has been accepted as an emerging subject which has its own standing based on analytical and experimental knowledge of the engineering disciplines concerned and numerical analysis. If we consider advances in fluid mechanical applications, we observe that major contributions have so far been rendered by a combination of experiments and approximate theoretical analysis. The analytical solution, albeit being a consequence of simplified form of governing equations, has the distinct advantage of immediately identifying some of the fundamental parameters of the given problem, and demonstrating some of the effects of these parameters on the physics of the problem. However, in order to include all the physical details of the problem formulation, total numerical simulation stepped in with its ability to handle the governing equations in their complete form during the end 1960's. Very soon it became a popular and reliable tool in engineering analysis. Today predictive procedures support experiments, enrich/extend the range of analytical solutions and finally contribute in product development.

Some of the major applications of numerical simulation are discerned in wind tunnel testing and combustion studies. In these applications the rapid decrease in the cost of computations (with the advent of modern high speed computers) compared to the increase in the cost of performing experiments with sophisticated gadgets has made numerical simulation an attractive alternative. The calculation of aerodynamic characteristics related to any new design through the application of numerical simulation is becoming indeed cheaper than measuring these characteristics in a wind tunnel. As such, in aircraft industries, the testing of preliminary design for new aircrafts or their components are first performed on the computer before taking-up rigorous wind tunnel tests. The wind tunnel is used to do the final fine-tuning of the design of new aircrafts and its components. In addition to economy, numerical simulation provides detailed flow field information when required. Even in very reliable experiments, measurements are generally taken at a few points. Over and above, many realistic conditions such as large sizes, very high temperature, toxic substances, fast transients which are indeed formidable to handle in experiments, can be simulated to some degree of confidence. Sometimes, in fluid mechanical studies, we are interested in investigating idealized conditions such as two dimensional or constant density flows. These tasks can be accomplished comfortably through numerical simulation.

Fig.1. Unstructured grid for two automobiles
(after Weatherill, 1988)


Fig. 2. Iso-surfaces of the spanwise vorticity (wz = + 0.25)
after Saha, 1999

An important question to address during the use of numerical simulations, however, is the validation of the predicted results. The results of numerical simulation are only as much as valid as the physical models incorporated in the governing equations. In addition, truncation errors associated with a particular algorithm employed to obtain a numerical are often amazingly accurate. Let us focus on some results generated by CFD. Figure 1 shows a typical computational domain and the grid-mesh used to discretize the domain for analysing flow around two automobiles (after Weatherill [47]). Figure 2 shows the isosurfaces of the spanwise vorticity in the wake of a square cylinder (courtesy Saha [48]).

It is also necessary to mention here than in many complex applications, experimental data base corroborate to a high degree accuracy, with numerical predictions. Keeping in mind such prediction power and acknowledging the fact that computations are frequently cheaper than conducting sophisticated experiments, engineers are more and more shifting towards the numerical simulations for testing and calibrating preliminary design.

Inspite of their extensive utility, prediction procedures are not completely devoid of limitations. For instance, they cannot reproduce physics that has not been properly included in the formulation of the problem. Thus, complex problems for which an adequate mathematical model support is not available can only be simulated through some approximate models. Accuracy of the results will definitely depend on the correctness of the model. For example, most of the solution procedures available for turbulent flow predictions use turbulence models which are indeed mere approximations of the real physical phenomena and the success of such models depend on the choice of some empirical constants. However, in recent times that there is also an emerging thought that on a fine enough scale, all turbulent flows obey the regular Navier Stokes equation and if a very fine grid can be used (order of 109 to 1012 grid points in a domain that typically has dimension of a few cms!) both the small scale and the large scale aspects of turbulence can be evaluated. This approach is known as Direct Numerical Simulation (DNS) of turbulence (Rai and Moin [49]).

Large Eddy Simulation (LES) is another technique intermediate between the direct simulation and the modeling approach discussed earlier. In LES, the contribution of the large scale structures to the momentum and energy transfer is computed exactly and the effect of the smallest scales of turbulence is modeled (Piomelli and Liu [50]). These developments indicate the close linkage between the level of physical understanding of a problem and the predictive procedures employed.

Another limitation of numerical simulations is discerned in reacting flows with moving flame fronts. Generally, there are a lot of uncertainties with respect to the determination of chemical reaction rates and the result of the numerical simulation of reactive flows are not free from such uncertainties. In addition, the presence of moving fronts complicate the problem tremendously. This is due to the fact that very small time steps are required to resolve the large gradient regions near the moving fronts. Since a numerical solution can only be as good as the grid employed to predict it, simulations of reactive flows result in the vicious circle of choosing appropriate grids which can adequately capture the important features of the unknown solutions. At times, if the choice of the numerical grid is not proper or the convergence of the iterative scheme is insufficient, the accuracy of the predictions may be unreliable. However in conclusion, it can be said that the results of numerical simulation are very much dependent on the physical understanding of a problem and the implementation of the subtle numerical considerations.

4. Closure

In the developed countries, the CFD has made rapid inroads in the academic and research institutions. Universities in the west and Japan possess, by definition, a synergy between teaching and research. They contribute to the society by two major products, professionals and technology. The CFD has brought about a major break-through in the technology-development in past twenty years. The Universities have been partly prodded into the technology development by generous support from the government agencies such as NASA (USA), DLR (Germany) and ONERA (France). In India, till date, the CFD is mostly practised in the IITs, IISc, DRDL (Hyderabad), VSSC (Trivandrum), NAL (Bangalore), ADA (Bangalore), ADE (Bangalore), GTRE (Bangalore), HAL (Bangalore), TRDDC (Pune) and BHEL (Hyderabad and Bhopal). Some well known multinational companies, such as FLUENT and ICEMCFD have started their operations in India. A need has been felt to promulgate the importance and utility of CFD to the engineering colleges (state and regional ) and medium-range industries. We hope that the next Millennium will benefit a low with such a spread of the CFD. The academia will come forward to help the industry in understanding the difficulties and the industry will expose the intricate problems of technology to the academia.


Van Dyke, M., Perturbation Methods in Fluid Mechanics, Academic Press (1964).

Howarth, L., Proceedings of Royal Society of London, A164, 547, 1938.

Schlichting, H., Boundary Layer Theory, McGraw Hill, New York, 1987.

Von Mises, R., Mathematical Theory of Compressible Fluid Flow, Academic Press, 1958.

Richardson, L.F., Phil. Trans. the Royal Society of London, Ser-A, 210, 307-357, 1910.

Liepman, H., Acad. Wiss., Math. Phys. Klasae, Sitzungsber, 3, 388, 1918.

Southwell, R.V., Relaxation Methods in Engineering Science, Oxford Univ. Press, 1940.

Frankel, S.P., Mathematical Tables and Other Aids to Computations, 4, 65-75, 1950.

Ritz, W., Zeitschrift fur Angewandte Mathematik and Mechanik, 35, 1-61, 1909.

Courant, R., Bulletin of the Americal Mathematical Society, 49, 1-23, 1943.

Clough, R.W., Proceedings of the Second ASCE Conference, Pittsburgh, 1961.

Courant, R., Friedrichs, K.O. and Lewy, H., Mathematische Annalen, 100, 32-74, 1928.

O'Brien, G.G., Hyman, M.A. and Kaplun, S., J. of Math. Phys., 29, 223-251, 1950.

Fromm, J.E., Los Alamos Scientific Laboratory Report, LA-2910, 1963.

Fromm, J.E., and Harlow, F.H., Phys. of Fluids, 6, 975-982, 1963.

Thoman, D.C. and Szewczyk, A.A., Phys. of Fluids, Supplement II, 76-87, 1969.

Rimon, Y. and Cheng, S.I., Phys. of Fluids, 12, 949-959, 1969.

Hamielec, A.E., Hoffman, T.W., and Ross, L.L., AIChE J, 13, 212-219, 1967.

Hamielec, A.E., Johnson, A.I. and Houghton, W.T., AIChE J, 13, 220-224, 1967.

Pearson, C.E., J. Fluid Mech., 21, 611-622, 1965.

Peaceman, D.W. and Rachford, H.H., J. Soc. Indust. Applied Math, 3, 28-41, 1955.

Douglas, J. and Rachford, H.H., Trans. Amer. Math. Soc., 82, 421-439, 1956.

Lax, P.D., Comm. Pure Appl. Math., 7, 159-193, 1954.

Evans, M.E. and Harlow, F.H., Los Alamos Report LA-2139, 1957.

Lax, P.D. and Wendroff, B., Comm. Pure Appl. Math., 13, 217-237, 1960.

MacCormack, R.W., AIAA Paper 69-354, Cincinnati, 1969.

Gary, J., Courant Institute Report, NY 09603, New York Univ, 1962.

Moretti, G. and Abbett, M., AIAA J. 4, 2136-2141, 1966.

Moretti, G. and Bleich, G., Sandia Report, SC-RR-68-3728, 1968.

Harlow, F.H. and Welch, J.E., Phys. Fluids, 8, 2182-2188, 1965.

Harlow, F.H. and Amsden, A.A., Los Alamos Report, LA-4370, 1970.

Chorin, A.J., Math. Comput., 22, 745-762, 1968.

Patankar, S.V., and Spalding, D.B., Int. J. Heat Mass Transfer, 15, 1787-1806, 1972.

Patankar, S.V., Numerical Heat Transfer and Fluid Flow, Hemisphere, 1980.

Van Doormaal, J.P. and Raithby, J.D., Numerical Heat Transfer, 7, 147-163, 1984.

Thompsan, J.F., Warsi, Z.U.A. and Mastin, C.W., Numerical Grid Generation, 1985.

Anderson, D.A., Advances in Grid Generation, ASME Fluids Eng. Conference, 1983.

Babuska, I., Chandra, J., and Flahertx, J.E. (Ed.), Adaptive Computational Methods for Partial Differential Equations, SIAM, 1983.

Baliga, B.R., PhD Thesis, Univ. of Minnesota, Mineeapolis, 1978.

Baliga, B.R. and Patnakar, S.V., Numerical Heat Transfer, 6, 245-261, 1983.

Baliga, B.R., Pham, T.T. and Patankar, S.V., Numerical Heat Transfer, 6, 263, 1983.

Taylor, C. and Hood, P., Computers and Fluids, 1, 73-100, 1973.

Jackson, C.P. and Cliffe, K.A., Int. J. Num. Meth. Engg., 17, 1980.

Taylor, C. and Hughes, T.G., Finite Element Programming of the Navier-Strokes Equations, Pineridge Press, 1987.

Peric, M., Ph.D Thesis, University of London, 1985.

Majumdar, S., Rodi, W. and Zhu, J., Journal of Fluids Engineering, 114, 496-503, 1992.

Weatherill, N.P., Int J Numerical Methods in Fluids, 8, 181-197, 1988.

Saha, A.K., Dynamical Characteristics of the Wake of a Square Cylinder at Low and High Reynolds Numbers, Ph.D. Thesis, Department of Mechanical Engineering, IIT Kanpur, 1999.

Rar, M.M. and Moin, P., J. Comput. Phys.,96, 15-53, 1991

Piomelli, V. and Liu, J., Phys. Fluids, 7, 839-848, 1995.

Gautam Biswas

Department of Mechanical Engineering

Indian Institute of Technology Kanpur

Kanpur- 208016

e.mail: gtm@iitk.ac.in

[back] [next]