CFD: Progress and Prospects

A lecture at the Asian Symposium ASCHT-2007 CFD: Progress and Prospects by Brian Spalding, of CHAM, Ltd

1. Introduction1.1 Purpose Computational fluid dynamics started half a century ago. In this lecture, I review its progress and seek to indicate how it may profitably develop further. I direct my words to research workersseeking problems which it is possible and beneficial to solve. I address also engineers, especially those working in process industries, whose designs can be improved if the indicated developments are carried out.

1.2 Patterns of analysis • The problems facing applied science are multi-dimensional; and they can be approached in various ways. • The main dimensions of variation are in: • time, • space, and • population (to be explained below). Variations in time are easiest to handle, because we all grow older at the same rate: one day per day. Variations in space are more complex, but easy to understand; for some of us can run faster than others.

1.2 Patterns of analysis Variations in population? Here is a one-dimensional histogram representing the distribution of the age of persons for a particular community at a particular time; and here is a picture to show that histograms can be two-dimensional. • Populations which are relevant to CFD include those of: • liquid droplets with differing diameters; • solid particles with differing velocities; • gas ‘fragments’ with differing compositions, or temperatures; and • radiation fluxes with differing directions.

1.2 Patterns of analysis • I shall further distinguish the three main approaches to non-uniformity, whether in time, space or population dimensions,namely: • neglect, • presume, • which means in effect, guess, and • calculate; and I shall argue that, in respect of calculation, the methods which are used for spatial variations can be applied to populationvariations also.

1.2 Patterns of analysis • I shall not argue that 'neglect' is always bad, or that 'calculate' is always best. • Indeed, most successful approaches are hybrid; thus: • even the most extreme of the calculators neglect something; and • nearly all presume rather than calculate some non-uniformities. • What is necessary is to make wise decisions about • what to neglect, • what to presume, • what to calculate, and • when to do each.

1.3 The structure of the lecture • In part 2, I shall explain my 3-dimension ~ 3-approach classification; and I shall illustrate it by way of examples from science and engineering. • In part 3, I shall recommend that CFD specialists should provide: • heat-exchanger designers with software based on less presumption and more calculation; • chemical-reactor operators with prediction tools which calculate the distribution of fluid fragments in composition space; and • mechanical engineers with computer codes which calculate the flow of fluids and the stresses in solids simultaneously.

2. Examples of engineering analysis2.1 Piston engines; space-direction variations The steam engine For this example, the 'neglect' approach is quite satisfactory, because the variations of steam temperature and pressure with position in the space above the piston are small at any instant of time.

2. Examples of engineering analysis2.1 Piston engines; space-direction variations Internal-combustion engines Here the 'neglect' approach is not satisfactory,because flames spread slowly. The 'presume' approach is best, especially when flame speed or spray burning rates are based on experimental observations. The 'calculate' approach, i.e. conventional CFD, is often employed; with limited success. Why? Because it neglects 'population' aspects of: (1) turbulent combustion and (2) droplet vaporisation.

2. Examples of engineering analysis 2.2 Simpler turbulent flows . The plane turbulent mixing layer; non-uniformity in space I start with the simplest of all turbulent flows; the plane mixing layer. The task is to predict the angle of the wedge-shaped layer of turbulent fluid at the edge of a jet injected into fluid at rest.

Shape functions and weighting functions The 'neglect' approach is not applicable here; for non-uniformity is of the essence. • The presumed-profile approach involves: • Guess the shapes of the velocity and effective-viscosity profies, e.g. as sloping or horizontal straight lines • Multiply the differential equations by weighting functions. • Integrate across the layer analytically. • Deduce the angle by algebra. • Advantage: quick and easy. • Disadvantage: accuracy is uncertain.

The plane turbulent mixing layer; the Finite-Volume Method The ‘calculate’ approach (version of Patankar and myself, 1967): • presumes only that the velocity profile is a histogram, with unknown column heights; • uses weighting functions of 1, i.e.none at all; • integrates across each histogram interval; • deduces the unknowns numerically. This is now known as the 'finite-volume' method' (FVM),the general form of its equations being: value in the volume = sum for all faces of coefficient * value in neighbour volume + sum of additional sources wherein the coefficients express diffusion and convection.

Other steady-state turbulent jets, wakes, plumes and boundary layers The early days of CFD; a condensed history • The FVM was soon applied to these flows • which: • had already been extensively studied experimentally, and by presumed-profile methods; • are 'parabolic' (i.e. downstream events do not influence upstream ones); • therefore permitted solution by 'marching' methods' on memory-scarce computers; • allowed turbulence models to be tested; • gave us confidence to extend the FVM to recirculating, three-dimensional, unsteady, compressible and chemically-reacting flows

2.3 Steady flow around solid bodies immersed in fluid streams Streamlined objects Before CFD, • aircraft design was based mainly on a 'neglect' approach, in that the variations of stagnation pressure • were neglected. The aerodynamic forces on the aircraft were then computed by way of ideal-fluid theory. • The effects of viscosity, and indeed turbulence were expressed by the supposition that the 'displacement thickness' of thin boundary layers enveloping wings and fuselage made these, in effect, rather thicker than they truly were. • The presumption approach was used, however, to calculate the displacement-thickness distribution; so the whole method can be characterised as being 'hybrid'.

Current practice • Now that CFD exists, • the calculation' approach is adopted for the whole of the space occupied by the fluid; which allows also the small regions of 'separated flow’ to be simulated. • However, an accurate calculation of the frictional forces on the solid surface can be made only by the use a very fine grid in the boundary layer; • so, for economy, some element of profile-presumption is retained, by way of wall functions.

Flows around and inside buildings • Before CFD, flow prediction was based on experiments with small geometrically similar physical models; • but this was unreliable , because the similarity criteria of Reynolds (viscosity) and Froude (buoyancy) could not both be satisfied. • Neither the neglect nor presume approaches had anything to offer. Therefore, engineers concerned with heating, ventilating, air-conditioning and fire-protection of buildings were among the first to turn to CFD.

Flows around and inside buildings • CFD has satisfied their requirements; and • it is for widely used for simulating fires in car-parks and other buildings; • BUT, for phenomena such as the fire-ball, it needs to take account of variations in hot-gas-population space.

2.4 Chemical-engineering equipment Heat exchangers; non-uniformities in space No designer can 'neglect' the temperature variations in heat exchangers. Instead, most guess them as being similar to that calculated for idealised counter-flow systems. Since they know that the flow patterns must differ, they multiply their calculated heat-transfer rates by correction factors like those on the right. But these are stillguesses, none the less.

Heat exchangers; non-uniformities in space (end) These presumption practices derive from the pre-CFD age. However, it was shown more than thirty years ago (by Patankar and myself, as it happens), that the calculate approach is practicable and indeed easy. It is strange therefore that most heat exchangers today are still based on presumption rather than calculation. Therefore, in section 3.1, below, I shall be recommending a change of practice.

Stirred chemical reactors, showing variations in both space and population The process: Many chemicals products are created by pumping feedstock materials (A and B) into a reactor vessel, where they are stirred together by a paddle, in order to react chemically. The task is to predict how the rate of production of C from reactants A and B depends upon the power consumed by stirring and the rate when mixed in a test-tube, where: rate/(concA*concB) = k_tube .

Stirred chemical reactorsVariations of time-averaged concentration Before CFD, the 'neglect' approach had to be used for variations with position; and it was not bad; for, if the stirring is vigorous enough, the time-average values ofconcA and concB will indeed be almost uniform. But what about moderate stirring? The 'presume' approach is not usable in this case; for no guidance exists as to what profiles should be presumed. Nowadays, CFD is employed; but it is not enough; for, if R_ave / (concA_ave * concB_ave)= k_reactor , it is found experimentally is that k_reactor is much less than k_tube. Why is this?

Stirred chemical reactorsVariations in population space The answer: non-uniformity in population space, also called unmixedness,shown here -> At any point in the reactor, fluid fragments of many different concentrations can be found. To calculate their time-average values, one must know for what proportion of time each is present. That means that one needs a probability-density function, like this ---> Can one calculate it? Yes, as I shall explain later; and for each location and stirring rate too. From it can be deduced the C- production rate.

Furnaces and other combustors; more variations in space and population • General description • A coal-fired furnace is a special kind of chemical reactor; and the processes taking place in it present a severe challenge to computer simulation, because of the importance of: • chemical reactions (coal pyrolysis, volatilisation, combustion, NOX formation) • solid-fluid interaction (diffusion of oxygen to the surface); • thermal radiation; and • particle-wall impact.

Furnaces and other combustorsVariations in position and population Which approach should be used for space variations? Only the calculate approach has any hope of representing the distributions of temperature, velocity, and pressure throughout the volume; and it has indeed been used for many years. • And for population non-uniformity? • Of coal-particle size: often neglected but sometimes presumed to vary in accordance with the empirical formula of Rosin and Rammler; • of radiation angle: often neglected ( in conduction model) sometimespresumed (in six-flux model) , and less often calculated (discrete-ordinates formulation); • of radiation wavelength: nearly always neglected; • of gas concentrations : nearly always neglected. • To recommend calculate for all would be too ambitious.

2.5 Simpler non-uniformities in population: droplet-size Vaporization of fuel sprays (in Diesels or gas turbines) consisting of droplets of various diameters, D, which change size at a rate governed by: - dD/dT = const * (1/D) * ln(1+B) where B, the driving force for mass transfer, depends upon (e.g.) local temperatures and other gas properties. This shows that droplets diminish in size at different rates, the smaller ones disappearing the more rapidly. . The task is to calculate the overall rate of vaporization. This necessitates knowing the droplet-size distribution at each location and each time.

Vaporization of a spray; droplet-size population The usual three ways are: 1. Neglect variations, i.e. suppose that all the droplets at a single location in the spray have the same diameter. 2. Presume that the profile is constant (e. g.) of Rosin-Rammler form, which cannot be very accurate. 3. Calculate the ordinates of the histogram by way of a standard finite-volume equation, with the source term dD/dT above. Use calculate if droplet size is critical, as in fire extinction.

The turbulent diffusion flame; fuel-air-ratio population Experimentally-observed unmixedness Hottel, Weddell and Hawthorne drew attention in 1949 to the 'unmixedness' of the gases in a flame produced by a jet of fuel gas injected into air. They measured finite time-average concentrations of both fuel and oxygen at the same location. That could never be found in a laminar flame. The first CFD analyses It was not until 1971 that the first attempt to simulate this unmixedness numerically was made, on the basis of a very simple profile presumption.

The turbulent diffusion flame; presumed fuel-air-ratio population The guess was that, at a point where the time-average fuel-air ratio was F, say, the gases actually present there had the ratio F+ g for half the time, and F- gfor the other half. Standard CFD calculated F easily. For g, a new differential equations was invented, having sources guessed as being proportional to gradients ofF- and velocity. This approach, when appropriate empirical constants were introduced, allowed turbulent diffusion flames to be simulated.

Confined pre-mixed flame; reactedness population In the turbulent diffusion flame, fuel and air enter separately, and must be mixed before chemical reaction can occur, at a rate limited by the rate of that mixing. I now consider a flow in which the fuel and air are mixed before they enter, at uniform and constant velocity, a plane-walled duct in which is placed a bluff-body 'flame- holder'. A turbulent wedge-shaped flame spreads across the duct, as the sketch indicates; and the profile of longitudinal velocity is roughly as shown. What then limits its rate? A different kind of mixing: that between burned and unburned gases.

Confined pre-mixed flame;the near-constancy of its angle • When first investigated, this flame showed some puzzling features, namely that the wedge angle was almost independent of: • inlet velocity • fuel-air ratio; • inlet temperature; • pressure; and • inlet turbulence intensity. • But why? • H.S. Tsien, while at CalTech, explained the shape of the profile; but what governed its angle remained a mystery. • We learned only later • non-uniformity in space depends on • non-uniformity in population.

Confined pre-mixed flame; the first population presumption The guessed profile The first idea, embodied in the so-called eddy-break-up model , was that the gas population consisted of two components, namely: (1) fragments of wholly un-burned gas which were too cold to burn; and (2) fragments of hot wholly-burned gas which also could not burn because either all the fuel or all the oxygen had been consumed. The histogram representing the presumed population therefore consisted of two spikes; and their relative heights dictated what would be measured as the time-average temperature.

Confined pre-mixed flame;collision between burned and unburned gas fragments These two elements of the population were thought of as colliding with one another and thereby producing sub-fragments of intermediate temperature and composition. These latter, being sufficiently hot and also containing reactants, could burn; and did so very rapidly, thereby increasing the height of the right-hand spike. Their actual concentration was considered, implicitly, to be negligibly small. The rate of collision per unit volume was guessed as proportional to the rate of dissipation of turbulence energy. This explained why the flame angle remained almost unchanged when the inflow velocity was increased.

Confined pre-mixed flame;the next presumed reactedness profile The four-fluid model The EBU, published in 1970, became very popular; so much so that 25 years passed before the obvious next step was taken;: to increase the number of presumed components from 2 to 4 ! Collisions between fluids 1 and 3 created fluid 2, 2 and 4 created fluid 3, 1 and 4 created fluid 2 and also fluid 3. Reaction of fluid 3 created fluid 4 at a chemistry- controlled rate. Fluids: 1 2 3 4

Confined pre-mixed flame;applications of the four-fluid model The chemistry-controlled step (fluid 3 creates fluid 4) explained: why: 1. the flame angle remained nearly constant, and 2. the flame could be suddenly extinguished by a velocity increase. The four-fluid model was used successfully for simulating flame spread in a baffled duct and for oil-platform explosion simulation. It has been little used; but it was the first step towards calculating the reactedness population,

From four fluids to many: the multi-fluid model In conventional CFD, we divide space and time into as many intervals as we need. Why not do the same for the reactedness at each point? The height of each column can then be deduced from a Finite-Interval equation’ like this: height of interval= sum for all faces of coefficient * height of neighbour interval + sum of additional sources + sum for all other intervals of coefficient * height of other interval )

What the terms in the finite-interval equation represent • In: height of interval= sum for all faces of coefficient * • height of neighbour interval + • the coefficients express rates of convection and diffusion, as in the the finite-volume equations of conventional CFD. • But in: sum for all other intervals of coefficient * • height of other interval • the coefficients express the physical and chemical processes: • collision between members of the fluid population, and • chemical conversion of one member into another. • The finite-interval method is thus merely a natural extension of the finite-volume method; and its equations can be solved in the familiar successive-substitution manner. • The calculation of population distributions is easy.

How material is distributed after collision Here is a diagram from one of the earliest publications. It depicts one of the possible hypotheses, called 'Promiscuous Mendelian'. The 'colliders' are treated as 'mother' and 'father’; and the word 'promiscous' implies that any two members of the population may collide. The word Mendelian, a reference to Gregor Mendel, the Austrian "father of modern genetics", implies that the offspring may appear with equal probability in any interval between those of the parents.

A calculated probability-density function This hypothesis has been embodied in the PHOENICS computer code. Here is one reactedness histogram, computed with its aid. As in the the eddy-break-up guess, there are indeed spikes at zero and unity reactedness; but calculation has shown that the intervals in-between are alsopopulated. Such probability distributions can to be computed for each location in the flame. Then the desired reaction rate for the whole flame can be deduced.

Application to gas-turbine combustion A three-dimensional gaseous-fuel combustor I show here one sector of a simple combustor proposed by Professor Wu Chung-Hua in the early days of PHOENICS.

Smoke formation rate is influenced by turbulent fluctutions Much later, I used this combustor to show how one must not neglect fluctuations of fuel-air ratio when predicting smoke formation. I used a 10-fluid model, with fuel-air-ratio as the population-defining attribute. Each cell had its own computed histogram The differences, although small. are significant when CFD is being used to optimise the design.

Concluding remarks for Part 2 It has been shown that: • variations in population space should not be neglected especially when chemical reaction is involved; • they can be presumed; • but it is better to calculate them. Why are not crowds of researchers pouring into this scarce-explored territory? Perhaps because they are waiting for less-timid crowds to do so first.

3. Recommendations3.1 To heat-exchanger designers So far, I have been discussing general ideas. Now I wish to make three specific recommendations. Current practice I have already mentioned that heat exchangers are still designed in the basis of presumption. A shell-and-tube heat exchanger looking like this (tubes not shown) can be expected to have a rather complex flow in the shell.

or Yet the software used by designers presumes that the flow in the shell can be conceptualized thus, and described by very few parameters. But why presume when one can calculate, as was shown to bepossible by the 35-year-old publication in which this image appeared?

3.1 To heat exchanger designers The solution • The solution is: • do notattempt to calculate the flow pattern between the tubes in detail, because current computers are not large or fast enough to handle the necessary fine grids except for a few tubes at a time. • Instead, use the space-averaged approach, with empirically-based formulae for: •  heat-transfer coefficients per unit volume, and •  friction factors per unit volume, • as functions of local Reynolds and Prandtl numbers. • Then solve the finite-volumeequations for (space-averaged) velocity, pressure, temperature for the shell- and tube-side fluids, treating both as interpenetrating continua, as is easily possible.

3.1 To heat exchanger designers The solution (contd) I now show some (not new) results for (the central plane of symmetry of) a particular shell-and-tube heat exchanger. (a) The shell-side velocity vectors, when calculated, appear thus (b) The consequential shell-side temperatures, are not,aspresumed, a succession of vertical stripes; although the calculated tube-side temperatures are (very nearly).

3.1 To heat exchanger designers The solution (end) (c) The conventional heat-exchanger-design packages presume that the shell-side, tube-side and overallheat-transfer coefficients are uniform throughout; but calculation reveals that they are not, as the next pictures clearly demonstrate. Corresponding non-uniformities are exhibited by the calculated Reynolds- and Prandtl-number values, and the temperature-dependent fluid properties, from which the heat-transfer coefficients have been computed.

Recommendation number 1 • My conclusions are… • that the conventional presumptions are evidently incorrect; • that therefore software which is based on them will generate unsafe designs; and • that the calculate approach, using experimentally-based data for the space-averaged heat-transfer and friction coefficients, is the only sound basis for the design of heat-exchanger equipment. I declared at the beginning that I had something to say to engineers. This first recommendation is addressed to them: Demand that the suppliers of your heat-exchanger-design software build into it the calculate approach.

3.2 To stirred-reactor designers and operators The calculation required by my first recommendation concerned non-uniformities in space. There are therefore many CFD specialists who will know how to implement it. My second concerns non-uniformities in population; experts in these are harder to find. The task is to predict how stirring-rate influences the conversion rate of reactants A and B into C in reactors of the kind which I discussed in Part 2.

3.2 To stirred-reactor designers and operators (contd) An example I turn to a ten-year old work [Ref ], in order to emphasise that the idea is not new, merely neglected. It concerns, for simplicity, reactants for which the rate constant measured in a laboratory test tube (i.e. k_lab) is very large. The geometry, and the body-fitted-co-ordinate grid used in the CFD calculation, are shown below.

3.2 To stirred-reactor designers and operators (contd) • But what about the mixture-ratio population grid? • Two distinct cases were considered, namely that: • the materials from the entering streams of reactants A and B were fully mixed at each point in the reactor, which would correspond to presuming • that its pdf was the single spike shown on the following diagram, and that • the amount of product C was as indicated by its horizontal location; • alternatively, at each point there could be found varying amounts of 'fluids' (in the multi-fluid sense) having one of eleven distinct mixture ratios, so that its pdf could be that of the histogram also shown there.

CFD: Progress and Prospects