Direct numerical simulations of multiphase flows are a rapidly growing field. In addition to continuous development of new and better numerical techniques and more extensive studies of the problems discussed so far in this Book, researchers are increasingly looking at new and more complex physical problems. In this chapter we examine briefly a few such extensions. We do not attempt to give an exhaustive list of all new applications of methods based on the “one-fluid” formulation of the fluid equations, but we hope that this introduction to the literature will be useful for our readers.

Additional fields and surface physics

Broadly speaking, new physics consists of new field equations, new surface effects, or both. Adding a new field is the simplest extension, but new fields often add new time- and length-scales. Mass transfer in liquids, for example, can lead to boundary layers that are much thinner than those resulting from either heat transfer or fluid motion. Resolving these boundary layers may introduce much more stringent resolution requirements than those necessary for the same problem in the absence of mass transfer.

The simplest new physics is probably heat transfer, where an advection/diffusion equation is solved for the temperature, and many authors have already studied multifluid problems involving heat transfer. The main complication is that large variations in the thermal conductivity across an interface can require fine grids and often it is better to use the harmonic mean of the conductivities at grid points where it is not defined, in the same way as for large differences in the viscosity (see Section 3.4).

The equations governing multiphase flows, where a sharp interface separates immiscible fluids or phases, are presented in this chapter. We first derive the equations for flows without interfaces, in a relatively standard manner. Then we discuss the mathematical representation of a moving interface and the appropriate jump conditions needed to couple the equations across the interfaces. Finally, we introduce the so-called “one-fluid” approach, where the interface is introduced as a singular distribution in equations written for the whole flow field. The “one-fluid” form of the equations plays a fundamental rôle for the numerical methods discussed in the rest of the book.

General principles

The derivation of the governing equations is based on three general principles: the continuum hypothesis, the hypothesis of sharp interfaces, and the neglect of intermolecular forces. The assumption that fluids can be treated as a continuum is usually an excellent approximation. Real fluids are, of course, made of atoms or molecules. To understand the continuum hypothesis, consider the density or amount of mass per unit volume. If this amount were measured in a box of sufficiently small dimensions ℓ, it would be a wildly fluctuating quantity (see Batchelor (1970), for a detailed discussion). However, as the box side ℓ increases, the density becomes ever smoother, until it is well approximated by a smooth function ρ. For liquids in ambient conditions this happens for ℓ above a few tens of nanometers (1 nm = 10−9 m).

Gas–liquid multiphase flows play an essential role in the workings of Nature and the enterprises of mankind. Our everyday encounter with liquids is nearly always at a free surface, such as when drinking, washing, rinsing, and cooking. Similarly, such flows are in abundance in industrial applications: heat transfer by boiling is the preferred mode in both conventional and nuclear power plants, and bubble driven circulation systems are used in metal processing operations such as steel making, ladle metallurgy, and the secondary refining of aluminum and copper. A significant fraction of the energy needs of mankind is met by burning liquid fuel, and a liquid must evaporate before it burns. In almost all cases the liquid is therefore atomized to generate a large number of small droplets and, hence, a large surface area. Indeed, except for drag (including pressure drops in pipes) and mixing of gaseous fuels, we would not be far off to assert that nearly all industrial applications of fluids involve a multiphase flow of one sort or another. Sometimes, one of the phases is a solid, such as in slurries and fluidized beds, but in a large number of applications one phase is a liquid and the other is a gas. Of natural gas–liquid multiphase flows, rain is perhaps the experience that first comes to mind, but bubbles and droplets play a major role in the exchange of heat and mass between the oceans and the atmosphere and in volcanic explosions.

The one-field formulation of the Navier–Stokes equations described in Chapter 2, where a single set of equations is used to describe the motion of all the fluids present, allows us to use numerical methods developed for single-phase flows. There are, however, two complications: the material properties (usually density and viscosity) generally vary from one fluid to the other and to set these properties we must construct an indicator function that identifies each fluid. We must usually also find the surface tension at the interface. The advection of the indicator function is the topic of Chapters 4 to 6 and finding the surface tension will be dealt with in Chapter 7. In this chapter we discuss numerical methods to solve the Navier–Stokes equations, allowing for variable density and viscosity. We will use the finite-volume method and limit the presentation to regular Cartesian grids. Since the multiphase flows considered in this book all involve relatively low velocities, we will assume incompressible flows.

For any numerical solution of the time-dependent Navier–Stokes equations it is necessary to decide:

1. (i) how the grid points, where the various discrete approximations are stored, are arranged;

2. (ii) how the velocity field is integrated in time;

3. (iii) how the advection and the viscous terms are discretized;

4. (iv) how the pressure equation, resulting from the incompressibility condition, is solved; and

5. (v) how boundary conditions are implemented.

These tasks can be accomplished in a variety of ways, but the approach outlined here has been widely used for multiphase flow simulations and results in a reasonably accurate and robust numerical method.

Progress is usually a sequence of events where advances in one field open up new opportunities in another, which in turn makes it possible to push yet another field forward, and so on. Thus, the development of fast and powerful computers has led to the development of new numerical methods for direct numerical simulations (DNS) of multiphase flows that have produced detailed studies and improved knowledge of multiphase flows. While the origin of DNS of multiphase flows goes back to the beginning of computational fluid dynamics in the early sixties, it is only in the last decade and a half that the field has taken off. We, the authors of this book, have had the privilege of being among the pioneers in the development of these methods and among the first researchers to apply DNS to study relatively complex multiphase flows. We have also had the opportunity to follow the progress of others closely, as participants in numerous meetings, as visitors to many laboratories, and as editors of scientific journals such as the Journal of Computational Physics and the International Journal of Multiphase Flows. To us, the state of the art can be summarized by two observations:

• Even though there are superficial differences between the various approaches being pursued for DNS of multiphase flows, the similarities and commonalities of the approaches are considerably greater than the differences. […]

Introduction

Understanding and predicting bubbly flows is of critical importance in a large number of industrial applications, including boiling heat transfer in power plants, various metallurgical processes, and in bubble columns in the chemical industry. In bubble columns, used for partial oxidation of ethylene to acetaldehyde, isobutene separation, wet oxidation of heavily polluted effluent, and the production of synthetic fuels, for example, gas is injected at the bottom and, as the bubbles rise, the gas diffuses into the liquid and reacts (Furusaki et al., 2001). Bubble columns ranging from tens to hundreds of cubic meters are common in the chemical industry and up to thousands of cubic meters in bioreactors, where longer process times are needed. The absence of any moving parts and their relatively simple construction makes bubble columns particularly attractive for large-scale operations (Deckwer, 1992). Their operation, however, is usually dependent on the size of the vessel and the difficulty of scaling up small pilot models makes numerical predictions important. Similar considerations apply to other bubble systems.

Computational modeling of industrial-size multiphase flow systems must by necessity rely on models of the average flow. Such models range from simple mixture models to more sophisticated two-fluid models, where separate equations are solved for the dispersed and the continuous phase. Since no attempt is made to resolve the unsteady motion of individual bubbles, closure relations are necessary for the unresolved motion and the forces between the bubbles and the liquid.

Droplet collisions and impacts are so spectacular that they have come to symbolize the beauty and fascination of fluid mechanics. Although simulations of two dimensional and axisymmetric systems go back to the early times of two-phase flow simulation, those of fully three-dimensional configurations have become possible only recently. It remains difficult, however, to perform realistic simulations of laboratory experiments.

Introduction

Droplet impacts are of major industrial interest. In what is perhaps the most significant application, fuel droplets impact on the walls of pipes and combustion chambers. There they may spread and form thin films or shatter into a spray of smaller droplets. Impacts also have an obvious relevance to ink-jet printing and spray coating. In other industrial processes, droplet impacts are of interest in metallurgy (Liow et al., 1996; Bierbrauer, 1995) and gas-injection processes. High speed droplet impacts may damage turbines operating with multiphase flows. In hypothetical severe nuclear reactor accidents, molten-core debris may impact on containment walls, splashing at very large velocity. In agriculture, impacts are related to the effect of rain on soil erosion (Farmer, 1973), or the spread of pesticides as they are sprayed on plants. Rain also influences air–sea interactions, enhancing the gas exchange and perhaps damping sea waves (Sainsbury and Cheeseman, 1950; Tsimplis and Thorpe, 1989). Droplet breakup, atomization, impacts, and splashes also cause the accumulation of charge in droplets, as shown by the 1905 Nobel physics laureate Philip Lenard following the work of Hertz (Lenard, 1892).

In Chapter 4 we saw several families of methods for locating and advecting the interface. Here, we focus on one of them: the VOF method. Historically, it was used for free surface flows with only one “fluid,” although it is now routinely used for two-fluid flows. In the VOF method the marker function is represented by the fraction of a computational grid cell which is occupied by the fluid assumed to be the reference phase.

A very large number (probably dozens) of VOF methods have been proposed. When we choose the method we try to strike a balance between several qualities: accuracy, simplicity and volume conservation.

Basic properties

The volume fraction or color function C is the discrete version of the characteristic function H; see Equation (4.3). We will be considering only two-phase or free-surface flows, so that the C data represent the fraction of each grid cell occupied by the reference phase. Furthermore, we restrict our analysis to Cartesian grids with square cells of side h = Δx = Δy.

The function C varies between the constant value one in full cells to zero in empty cells, while mixed cells with an intermediate value of C define the transition region where the interface is localized.

Low-order VOF methods do not need to specify the location of the interface in the transition region, but a geometrical interpretation of these methods shows that in two dimensions the interface line in each mixed cell is represented by a segment parallel to one of the two coordinate axes.

When the governing equations are solved on a fixed grid, using one set of equations for the whole flow field, the different fluids must be identified in some way. This is generally done by using a marker function that takes different values in the different fluids. Sometimes a material property, such as the fluid density for incompressible fluids, can serve as a marker function, but here we shall assume that the rôle of the marker function is only to identify the different fluids. As the fluids move, and the boundary between the different fluids changes location, the marker function must be updated. Updating the marker function accurately is both critical for the success of simulations of multiphase flows and also surprisingly difficult. In this chapter we discuss the difficulties with advecting the marker function directly and the various methods that have been developed to overcome these difficulties.

The VOF method is the oldest method to advect a marker function and – after many improvements and innovations – continues to be widely used. Other marker function methods include the level-set method, the phase-field method, and the CIP method. Instead of advecting the marker function directly, the boundary between the different fluids can also be tracked using marker points, and the marker function then reconstructed from the location of the interface. Methods using marker points are generally referred to as “front-tracking” methods to distinguish them from “front-capturing” methods, where the marker function is advected directly.

Various applications and natural processes involve large deformations and eventual breakup of liquid jets, layers, and droplets. When liquid masses fragment in a small number of pieces one speaks of breakup. More intense phenomena where, for instance, a liquid jet is broken into seemingly microscopic droplets are called atomization, although the term is somewhat incorrect, since the individual pieces are still far larger than atomic scales.

Nevertheless, atomization is a striking process in which finely divided sprays or droplet clouds are produced. This is often based on the ejection of a high-speed liquid jet from an atomizer nozzle. Many other configurations exist, such as sheets ejected at high speed from diversely shaped nozzles, or colliding with each other. As with many of the multiphase phenomena investigated in this book, atomization offers a rich physical phenomenology which is still poorly understood. Considerable progress has been made in the development of methods for atomization simulations during the last few years, and advances in hardware are making it possible to conduct simulations of unprecedented complexity.

Introduction

There are many important motivations for the study of spray formation, droplet breakup, and atomization. To take a first example from natural phenomena, spray formation atop ocean waves occurs when sufficiently strong winds strip droplets from the crests of the waves. Breaking waves also create bubbles that, when bursting at the surface, create a very fine mist that can rise high into the atmosphere.

The accurate computation of the surface tension is perhaps one of the most critical aspects of any method designed to follow the motion of the boundary between immiscible fluids for a long time. In methods based on the “one-fluid” formulation, where a fixed grid is used to find the motion, surface tension is added as a body force to the discrete version of the Navier–Stokes equations. Finding the surface force depends on whether the fluid interface is tracked by a direct advection of a marker function (VOF) or whether discrete points are used to mark the interface (front tracking). Here, we describe how to find surface tension for both VOF and front-tracking methods. At the end of the chapter we examine the performance of the various methods and the challenges in computing surface tension accurately and robustly.

Computing surface tension from marker functions

A marker function such as the color function C of the VOF method or the level-set function F is a function that indicates (marks) where the interface is. A whole family of methods for surface tension have been developed for the use of marker functions; however, these methods may also be used in connection with front tracking methods. The standard approach is termed the continuous surface force (CSF) method. We also describe a variant that conserves momentum exactly, the continuous surface stress (CSS) method.

Continuous surface force method

For simplicity, we consider first the case where σ is constant.

Instead of advecting a marker function identifying the different fluids directly, as in the VOF method, the boundary between the fluids can be represented by connected marker points that are moved by the fluid. This approach is usually called front tracking, and in Chapter 4 we discussed the basic idea briefly and gave a short historical overview. In this chapter we describe front tracking in more detail.

The use of connected marker points to track the motion of a complex and deforming fluid interface can lead to several different methods, depending on the details of the implementation. Generally, however, front tracking involves the following considerations:

1. (i) The data structure used to describe the front. Although the use of marker particles simplifies many aspects of the advection of a fluid interface, other aspects become more complex. We use front to refer to the complete set of computational objects used to represent the interface. In addition to the marker points, the front often includes information about the connectivity of the points, as well as a description of the physics at the interface. The management of the front can be greatly simplified by the use of the appropriate data structure. There is a fundamental difference in the level of complexity between fronts in two dimensions and in three dimensions and, in general, any data structure can be made to work reasonably efficiently in two dimensions. […]