Historical Perspectives of CFD and Some Applications

Dr. Gautam Biswas
Department of Mechanical Engineering, Indian Institute of Technology Kanpur and
Director, Central Mechanical Engineering Research Institute (CSIR), Durgapur-713209

1. Introduction

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 phenomena in space and 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 numbers. This is contrast to the symbolic expression of closed form analytical solution. The final objective of most engineering investigations is to obtain a quantitative description of the problem. 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 requires the repetitive manipulation of 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 and heat transfer problems.

2. Historical Perspective

The semi-analytical techniques which found wide use in fluid dynamics research were the perturbation methods, the similarity approach, and the integral methods, all for the viscous boundary layer calculations, and the methods of characteristics, 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. For structural problems involving elastic deformations, W. Ritz developed a method in 1909, 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. R. Courant in 1943 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 R.W. Clough in 1960. Since then, the method has made rapid strides for the modelling of structural engineering problems. The fluid flow and heat transfer modelling have been accomplished successfully in the recent years. A pioneering work on the uniqueness and existence of numerical solutions to partial differential equations was presented by R. Courant, K.O. Friedrichs and H. Lewy in 1928. The stability requirement (CFL condition) for hyperbolic partial differential equations was first shown in this work. The stability criteria for parabolic time-marching problems in Computational Fluid Dynamics were developed by Von Neumann. 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, F.H. Harlow developed an explicit forward time difference method at Los Alamos. Their method was used by D.C. Thoman and A.A. Szewczyk for cross flow over cylinders and Y. Rimon and S.I. Cheng for uniform flow over a sphere. The investigations were published in Physics of Fluids. Significant development was witnessed in the fifties and sixties towards the solution of inviscid compressible flow equations. Starting with shock-capturing technique of P.D. Lax 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. In 1960, a second order accurate finite difference scheme which reduces the shock smearing effect was proposed by P.D. Lax and B. Wendroff. This scheme later led to the development of the well known McCormack method. Although in the early simulation methods for viscous incompressible flow, vorticity and stream function were the calculated variables, in the sixties, simulations in terms of primitive variables (velocity components and pressure) began. Pioneering work in this direction was performed by F.H. Harlow and J.E. Welch and A.A Amsden at Los Alamos. These authors introduced explicit transient algorithms such as MAC and SMAC. Alexander Chorin 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 S.V. Patankar and D.B. Spalding at the imperial College. Based on this well-known SIMPLE algorithm and its later improvements such as SIMPLER and SIMPLEC, 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 discussions on these methods are provided in the book by J.F. Thompson , Z.U.A. Warsi and C.W. Mastin. Milovan Peric has proposed a general flow solver based on the application of conservation laws to non-orthogonal control volumes. Research is still being conducted on these methods and several complex flow simulations are being simulated using them.

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 cognate engineering disciplines 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. 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 and extend the range of analytical solutions and finally contribute in product development. Some of the major applications of numerical simulation have been discerned in wind tunnel testing and combustion studies. In these applications, the rapid decrease in the cost of computations 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. Over and above, many realistic conditions such as large sizes, very high temperatures, toxic substances, fast transients which are indeed formidable to handle in experiments, can be simulated to some degree to confidence. Let us focus on some results generated using CFD by the author and his co-workers.

3.1 Computation of Flow Structure and Heat Transfer due to Longitudinal Vortices in Heat Exchanger Applications

Longitudinal vortices have enormous utility for flow control. Longitudinal vortices are also capable of producing beneficial effects in transport enhancement. The vortices disrupt the growth of the boundary layer and serve ultimately to bring about enhancement of heat transfer between the fluid and its neighbouring surface. The present investigations determine the flow structure, in detail, behind the winglet type vortex generators placed in fully developed laminar and turbulent channel flows. Figure 3.1 shows temperature distribution on a fin surface in a fin-tube heat exchanger in the presence of a delta winglet pair. The complex flow structure in the flow passage formed by two neighbouring fins is revealed in Figure 3.2.

Fig. 3.1 Temperature distribution on a horizontal plane (at Y=0.0416 from the bottom plate)

Fig. 3.2 Cross-stream velocity at different axial locations from the inlet

Fig. 3.3 Vorticity contours on a cross plane behind a winglet (a) experiment (b) computation

The flow structure behind a winglet is complex and consists of a main vortex system, a corner (horseshoe) vortex system, and induced vortices. Hot-wire measurements were performed in order to corroborate the numerical predictions of the flow structure.

Fig. 3.4 Isolines for turbulent kinetic energy at different cross-sections for ReH/2 = 67,000

Figure 3.3 demonstrates a close comparison between the experimental and computation flow fields. Figure 3.4 compares the predicted kinetic energy contours with the experiments conducted by a well-known research group in Stanford University. The purpose of this study is to analyse performance of the delta winglet type vortex generators in improving heat transfer. Such vortex generators also show great promise for enhancing the heat transfer rate in plate-fin crossflow heat exchangers. The investigation was extended to compute flow and heat transfer in the air-cooled condensers for the Geothermal Power Plants.

Fig. 3.5 Nusselt number distribution on bottom plate of a channel with built-in circular tube and delta winglet-pair in common-flow-up configuration

This was a tripartite international project between Idaho National Laboratory (USA), Yokohama National University (Japan) and IIT Kanpur. The NEDO (Japan) was the sponsoring agency. The heat transfer distribution on the fin surface of a proposed heat exchanger module is shown in Figure 3.5. The solution was obtained following an algorithm developed by Prof. V. Eswaran. The winglets with common-flow-up configuration were found more beneficial in such applications.

3.2 Three-Dimensional Analysis of Flow in the Spiral Casing of a Reaction Turbine Using a Differently Weighted Petrov-Galerkin method

A weighted residual based finite element method was used to predict the flow structure and the head loss in the spiral casing of a reaction turbine. An explicit Eulerian velocity correction scheme was employed to solve the complete Reynolds-averaged Navier-Stokes equations. A Differently Weighted Petrov Galerkin (DWPG) method was used for spatial discretization. The simulation was performed for high Reynolds numbers. The head loss in the spiral casing was calculated for Re = 5 X 105 and Re = 5 X 106. The velocity field and the pressure distributions inside the spiral casing corroborate the trends available in literature.

Fig. 3.6 (a) Static pressure distribution on the horizontal mid-plane for Re = 106

A strong secondary flow is evolved on the cross-stream plane due to the centrifugal forces, which act at right angles to the main flow. Figures 3.6 (a) and (b) depict the pressure distribution and the secondary flow in the spiral casing of a hydraulic turbine.

Fig. 3.6 (b) Secondary flow at different cross sections for Re = 106 (a) θ = 0o; (b) θ = 90o; (c) θ = 180o; (d) θ = 240o

3.3 Flow Past Bluff Bodies

The flow past bluff body involves very complex flow physics as there are three simultaneous instabilities present in such flows, namely boundary layer instability, separated shear layer instability and primary vortex or Karman vortex instability. These three different frequencies interact among themselves and produce a complex flow structure. At the transitional Reynolds number (transition to three dimensionality), two distinct modes of vortex shedding are observed which are characterized by the spanwise wavelengths of the streamwise vortices.

The early modes at the transitional Reynolds number, called Mode-A, has a wavelength of three-cylinder width and Mode-B which occurs at the later part of the transitional regime has a wavelength of one-cylinder width. Mode-A also shows an intermittent low frequency irregularity in the shedding cycle called Vortex Dislocation. Figure 3.7 captures the temporal evolution of vortex dislocation in the wake of a square cylinder.

Fig. 3.7 Temporal evolution of vortex dislocation at various time t: (a) 510.11, (b) 541.99, (c) 553.81, (d) 561.59 and (e) 572.23 for a Reynolds number of 175

3.4 Large-eddy Simulation of Flow and Heat Transfer in Impinging Slot Jets

The flow field (Figure 3.8) due to an impinging jet at a moderately high Reynolds number, emanating from a rectangular slot nozzle was computed using Large Eddy Simulation (LES) technique. A dynamic subgrid-scale model was used for the small scales of turbulence.

Fig. 3.8 Instantaneous velocity field on x-z plane, ReB = 5800

Quite a few successful applications of the dynamic subgrid-scale stress model use planer averaging to avoid ill conditioning of the model coefficient. However, a novel localization procedure has been attempted herein to find out the spatially varying coefficient of the flow. The flow field is characterized by entrainment at the boundaries. Periodic boundary conditions could not be used on all the boundaries. The results reveal the nuances of the vortical structures that are characteristic of jet flows. The stress budget also captures locally negative turbulence production rate. Figure 3.9 displays a periodic variation of a temperature field due to radial jet reattachment.

Fig. 3.9 Thermal field for a radial jet reattachment flow

3.5 Simulation of Boiling Flows

Figure 3.10 shows simulation of film boiling and bubble formation in water at near critical conditions on an isothermal horizontal surface using a Coupled Level Set and Volume of Fluid (CLSVOF) based interface tracking method.

Fig. 3.10 Bubble release from the constant wall temperature surface

The complete Navier-Stokes equations and thermal energy equations were solved in conjunction with an interface mass transfer model.

4. Closure

A numerical solution can only be as good as the grid employed to predict it. Generation of appropriate grids, which can adequately capture the important features of the solutions, is a very important aspect. The results of numerical simulation are very much dependent on the physical understanding of a problem and the implementation of the subtle numerical considerations. A common numerical difficulty of applying CFD to flows at high Reynolds numbers and on rough walls is the resolution of the velocity gradients within the layer of diminishing thickness. This has prompted to use wall functions, which are often criticised for introducing errors in complex flows. Low-Reynolds number models attempt to reproduce the experimentally observed viscous layer scaling laws, based on molecular viscosity and turbulent friction velocity. However, one should not expect them to perform well near separation. The Large Eddy Simulation (LES) technique has demonstrated superiority over the RANS approach. In many LES simulations, it was possible to go very near to the wall (y+~1). LES requires excessive CPU time and memory for computations of 3D complex flow problems particularly at large Reynolds numbers. Numerical experiments based on LES and attempts to determine the quantities that are difficult to measure in Laboratories will have a critical role in future research in Computational Fluid Dynamics.


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  36. Thompson, J.F., Warsi, Z.U.A. and Mastin, C.W., Numerical Grid Generation (1985).

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

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

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

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

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

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

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

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

  45. Peric, M., Ph.D Thesis, University of London (1985).

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

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

  48. Saha, A. K., Biswas, G., and Muralidhar, K., Int. J. Heat and Fluid Flow, Vol. 24, pp. 54-66 (2003).

  49. Rai, M.M. and Moin, P., J. Comput. Phys.,96, 15-53 (1991).

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

  51. Biswas, G., Mitra, N. K., and Fiebig, M., Int. J. of Heat and Mass Transfer, vol.37, pp. 283-291 (1994).

  52. Maji, P. K. and Biswas, G., Int. J. for Numerical Methods in Engineering, Vol. 45, pp. 147-174 (1999).

  53. Cziesla, T., Biswas, G., Chattopadhyay, H., and Mitra, N. K., Int. J. Heat and Fluid Flow, Vol. 22, pp. 500-508 (2001).

  54. Jain, A., Biswas, G., and Maurya, D., Numerical Heat Transfer, Part A, Vol. 43, pp. 201-219 (2003).

  55. Tomar, G., Biswas, G., Sharma, A., and Agrawal, A., Phys. of Fluids, Vol. 17, 112103 (2005).

  56. Gerlach, D., Tomar, G., Biswas, G., and Durst, F., Int. J. of Heat and Mass Transfer, Vol. 49, pp. 740-754 (2006).

  57. Tomar, G., Biswas, G., Sharma, A., and Welch, S. W. J., Phys. of Fluids, Vol. 21, 032107-1 - 032107-8 (2009).

  58. Ray, B., Biswas, G., and Sharma, A., J. of Fluid Mechanics, Vol. 655, pp. 72-104 (2010).