A liquid foam is a dispersion of gas bubbles in a liquid matrix containing surface-active agents. Its flow involves the relative motion of bubbles, which switch neighbours during a so-called topological rearrangement of type 1 (T1). The dynamics of T1 events, as well as foam rheology, have been extensively studied, and experimental results point to the key role played by surfactants in these processes. However, the complex and multiscale nature of the system has so far impeded a complete understanding of the mechanisms involved. In this work, we investigate numerically the effect of surfactants on the rheological response of a 2D sheared bubble cluster. To do so, a level-set method previously employed for simulation of two-phase flow has been extended to include the effects of surfactants. The dynamical processes of the surfactants – diffusion in the liquid and along the interface, adsorption/desorption at the interface – and their coupling with the flow – surfactant advection and Laplace and Marangoni stresses at the interface – are all taken into account explicitly. Through a systematic study of the Biot, capillary and Péclet numbers that characterise the surfactant properties in the simulation, we find that the presence of surfactants can affect the liquid/gas hydrodynamic boundary condition (from a rigid-like situation to a mobile one), which modifies the nature of the flow in the volume from a purely extensional situation to a shear. Furthermore, the work done by surface tension (the 2D analogue of the work by pressure forces), resulting from surfactant and interface dynamics, can be interpreted as an effective dissipation, which reaches a maximum for a Péclet number of order unity. Our results, obtained at high liquid fraction, should provide a reference point, with which experiments and models of T1 dynamics and foam rheology can be compared.

We investigate the dispersion of a passive scalar such as the concentration of a chemical species, or temperature, in homogeneous bubbly suspensions, by determining an effective diffusivity tensor. Defining the longitudinal and transverse components of this tensor with respect to the direction of averaged bubble rise velocity in a zero mixture velocity frame of reference, we focus on the convective contribution thereof, this being expected to be dominant in commonly encountered bubbly flows. We first extend the theory of Koch et al. ( J. Fluid Mech. , vol. 200, 1989, pp. 173–188) (which is for dispersion in fixed beds of solid particles under Stokes flow) to account for weak inertial effects in the case of ordered suspensions. In the limits of low and of high Péclet number, including the inertial effect of the flow does not affect the scaling of the effective diffusivity with respect to the Péclet number. These results are confirmed by direct numerical simulations performed in different flow regimes, for spherical or very deformed bubbles and from vanishingly small to moderate values of the Reynolds number. Scalar transport in arrays of freely rising bubbles is considered by us subsequently, using numerical simulations. In this case, the dispersion is found to be convectively enhanced at low Péclet number, like in ordered arrays. At high Péclet number, the Taylor dispersion scaling obtained for ordered configurations is replaced by one characterizing a purely mechanical dispersion, as in random media, even if the level of disorder is very low.

We study unstable waves in gas–liquid two-layer channel flows driven by a pressure gradient, under stable stratification, not assumed to be set in motion impulsively. The basis of the study is direct numerical simulation (DNS) of the two-phase Navier–Stokes equations in two and three dimensions for moderately large Reynolds numbers, accompanied by a theoretical description of the dynamics in the linear regime (Orr–Sommerfeld–Squire equations). The results are compared and contrasted across a range of density ratios $r=\unicode[STIX]{x1D70C}_{liquid}/\unicode[STIX]{x1D70C}_{gas}$. Linear theory indicates that the growth rate of small-amplitude interfacial disturbances generally decreases with increasing $r$; at the same time, the cutoff wavenumbers in both streamwise and spanwise directions increase, leading to an ever-increasing range of unstable wavenumbers, albeit with diminished growth rates. The analysis also demonstrates that the most dangerous mode is two-dimensional in all cases considered. The results of a comparison between the DNS and linear theory demonstrate a consistency between the two approaches: as such, the route to a three-dimensional flow pattern is direct in these cases, i.e. through the strong influence of the linear instability. We also characterize the nonlinear behaviour of the system, and we establish that the disturbance vorticity field in two-dimensional systems is consistent with a mechanism proposed previously by Hinch (J. Fluid Mech., vol. 144, 1984, p. 463) for weakly inertial flows. A flow-pattern map constructed from two-dimensional numerical simulations is used to describe the various flow regimes observed as a function of density ratio, Reynolds number and Weber number. Corresponding simulations in three dimensions confirm that the flow-pattern map can be used to infer the fate of the interface there also, and show strong three-dimensionality in cases that exhibit violent behaviour in two dimensions, or otherwise the development of behaviour that is nearly two-dimensional behaviour possibly with the formation of a capillary ridge. The three-dimensional vorticity field is also analysed, thereby demonstrating how streamwise vorticity arises from the growth of otherwise two-dimensional modes.

A computational method for the simulation of viscous and compressible gas–gas flows is presented. It consists in solving the Navier–Stokes equations associated with a convection equation governing the motion of the interface between two gases using high-order finite-difference schemes. A discontinuity-capturing methodology based on sensors and a spatial filter enables capturing shock waves and deformable interfaces. One-dimensional test cases are performed as validation and to justify choices in the numerical method. The results compare well with analytical solutions. Shock waves and interfaces are accurately propagated, and remain sharp. Subsequently, two-dimensional flows are considered including viscosity and thermal conductivity. In Richtmyer–Meshkov instability, generated on an air–SF6 interface, the influence of the mesh refinement on the instability shape is studied, and the temporal variations of the instability amplitude is compared with experimental data. Finally, for a plane shock wave propagating in air and impacting a cylindrical bubble filled with helium or R22, numerical Schlieren pictures obtained using different grid refinements are found to compare well with experimental shadow-photographs. The mass conservation is verified from the temporal variations of the mass of the bubble. The mean velocities of pressure waves and bubble interface are similar to those obtained experimentally.

A commonly used class of methods for the numerical simulation of two-phase flows is level set. It is often reported though that this method does not accurately conserve mass of each fluid, unlike other interface capturing techniques such as volume-of-fluid. A further concern besides mass conservation is the formation of any parasitic currents. Since the initial formulation of level-set methods, however, numerous modifications have been proposed, and it does not seem clear whether mass conservation errors and parasitic currents are problematic for all of these and, if not, what key steps could be taken to avoid them. Furthermore, results reported in the literature are often for benchmark tests in two dimensions, and it is not clear whether a good performance there holds up in three dimensions. We undertake here a comparative study, reporting test results in two and three dimensions for various level-set methods on a variety of problems. Kinematical tests are first performed for prescribed velocity fields, followed by benchmark tests including the solution of the Navier–Stokes equations. It is shown that higher-order schemes for spatial and temporal discretization may improve mass conservation and avoid interface distortion. In particular, two reinitialization methods that are straightforward to implement perform very well at all these tests. It is demonstrated that some schemes introduce parasitic oscillations in the simulation of Rayleigh–Taylor instability.

Various expressions have been proposed previously for the rise velocity of gas bubbles in homogeneous steady bubbly flows, generally a monotonically decreasing function of the bubble volume fraction. For suspensions of freely moving bubbles, some of these are of the form expected for ordered arrays of bubbles, and vice versa, as they do not reduce to the behaviour expected theoretically in the dilute limit. The microstructure of weakly inhomogeneous bubbly flows not being known generally, the effect of microstructure is an important consideration. We revisit this problem here for bubbly flows at small to moderate Reynolds number values for deformable bubbles, using direct numerical simulation and analysis. For ordered suspensions, the rise velocity is demonstrated not to be monotonically decreasing with volume fraction due to cooperative wake interactions. The fore-and-aft asymmetry of an isolated ellipsoidal bubble is reversed upon increasing the volume fraction, and the bubble aspect ratio approaches unity. Recent work on rising bubble pairs is used to explain most of these results; the present work therefore forms a platform of extending the former to suspensions of many bubbles. We adopt this new strategy also to support the existence of the oblique rise of ordered suspensions, the possibility of which is also demonstrated analytically. Finally, we demonstrate that most of the trends observed in ordered systems also appear in freely evolving suspensions. These similarities are supported by prior experimental measurements and attributed to the fact that free bubbles keep the same neighbours for extended periods of time.

A new method is presented to perform three-dimensional simulations of two-phase flows with moving contact lines using level-set. To account for the full range of length scales involved in the physical problem under realistic conditions, without having to resolve the flow down to the smallest continuum scale, a dynamic contact angle model based on asymptotic theory is used in conjunction with the computational method. Contact-angle hysteresis is also represented in this methodology. The method is validated against simulations wherein the flow is fully resolved over all length scales, and experiments of spreading droplets and droplets sliding down an inclined substrate.

Viscous simulations of shock-bubble interaction and Richtmyer-Meshkov instability are performed using an explicit high-order computational method. The simulations are performed by solving the NavierStokes equations associated with two convection equations governing the interface between two fluids. The stiffened equation of state is used to relate the pressure to the total energy of a liquid or a gas. Two-dimensional two-phase flows are considered. The first flow concerns the Richtmyer-Meshkov instability developed on a post-shocked interface between air and sulphur hexafluoride (SF6). The influence of the grid refinement on the instability shape is studied. The second problem deals with a shock wave propagating in air and hitting a cylindrical bubble filled with helium or chlorodifluoromethane (R22). A spatial-time diagram represents the locations of the various pressure waves generated from the interaction between the shock wave and the interface. For both simulations, the numerical results are in agreement with experimental data and visualizations.

Axisymmetric non-isothermal spreading/dewetting of droplets on a substrate is studied, wherein the surface tension is a function of temperature, resulting in Marangoni stresses. A lubrication theory is first extended to determine the drop shape for spreading/dewetting limited by slip. It is demonstrated that an apparent angle inferred from a fitted spherical cap shape does not relate to the contact-line speed as it would under isothermal conditions. Also, a power law for the thermocapillary spreading rate versus time is derived. Results obtained with direct numerical simulations (DNS), using a slip length down to $O(10^{-4})$ times the drop diameter, confirm predictions from lubrication theory. The DNS results further show that the behaviour predicted by the lubrication theory – that a cold wall promotes spreading, and a hot wall promotes dewetting – is reversed at sufficiently large contact angles and/or viscosity of the surrounding fluid. This behaviour is summarized in a phase diagram, and a simple model that supports this finding is presented. Although the key results are found to be robust when accounting for heat conduction in the substrate, a critical thickness of the substrate is identified above which wall conduction significantly modifies wetting behaviour.

Computational methods have been extended recently to allow for the presence of moving contact lines in simulated two-phase flows. The predictive capability offered by these methods is potentially large, joining theoretical and experimental methods. Several challenges rather unique to this area need to be overcome, however, notably regarding the conditions near a moving contact line and the very large separation of length scales in these flows. We first summarize the main models for moving contact lines and follow with an overview of computational methods that includes direct continuum approaches and macroscale models that resolve only the large-scale flow by modeling the effects of the conditions near the contact line using theory. Results are presented for contact-line motion on ideal as well as patterned and grooved surfaces and for extensions to account for complexities such as thermocapillarity and phase change.

Abstract We consider the linear and nonlinear stability of two-phase density-matched but viscosity-contrasted fluids subject to laminar Poiseuille flow in a channel, paying particular attention to the formation of three-dimensional waves. A combination of Orr–Sommerfeld–Squire analysis (both modal and non-modal) with direct numerical simulation of the three-dimensional two-phase Navier–Stokes equations is used. For the parameter regimes under consideration, under linear theory, the most unstable waves are two-dimensional. Nevertheless, we demonstrate several mechanisms whereby three-dimensional waves enter the system, and dominate at late time. There exists a direct route, whereby three-dimensional waves are amplified by the standard linear mechanism; for certain parameter classes, such waves grow at a rate less than but comparable to that of the most dangerous two-dimensional mode. Additionally, there is a weakly nonlinear route, whereby a purely spanwise wave grows according to transient linear theory and subsequently couples to a streamwise mode in weakly nonlinear fashion. Consideration is also given to the ultimate state of these waves: persistent three-dimensional nonlinear waves are stretched and distorted by the base flow, thereby producing regimes of ligaments, ‘sheets’ or ‘interfacial turbulence’. Depending on the parameter regime, these regimes are observed either in isolation, or acting together.

Computational methods have been extended recently to allow for the presence of moving contact lines in simulated two-phase flows. The predictive capability offered by these methods is potentially large, joining theoretical and experimental methods. Several challenges rather unique to this area need to be overcome, however, notably regarding the conditions near a moving contact line and the very large separation of length scales in these flows. We first summarize the main models for moving contact lines and follow with an overview of computational methods that includes direct continuum approaches and macroscale models that resolve only the large-scale flow by modeling the effects of the conditions near the contact line using theory. Results are presented for contact-line motion on ideal as well as patterned and grooved surfaces and for extensions to account for complexities such as thermocapillarity and phase change.

AbstractWe study two-phase stratified flow where the bottom layer is a thin laminar liquid and the upper layer is a fully developed gas flow. The gas flow can be laminar or turbulent. To determine the boundary between convective and absolute instability, we use Orr–Sommerfeld stability theory, and a combination of linear modal analysis and ray analysis. For turbulent gas flow, and for the density ratio$r= 1000$, we find large regions of parameter space that produce absolute instability. These parameter regimes involve viscosity ratios of direct relevance to oil and gas flows. If, instead, the gas layer is laminar, absolute instability persists for the density ratio$r= 1000$, although the convective/absolute stability boundary occurs at a viscosity ratio that is an order of magnitude smaller than in the turbulent case. Two further unstable temporal modes exist in both the laminar and the turbulent cases, one of which can exclude absolute instability. We compare our results with an experimentally determined flow-regime map, and discuss the potential application of the present method to nonlinear analyses.

We propose an efficient level-set approach for numerical simulation of moving contact lines. The main purpose is to formulate and test a model wherein the macroscale flow is resolved while the effects of the microscopic region near a contact line are represented using asymptotic theories. The model covers viscous as well as inertial regimes. Test simulations include axisymmetric displacement flow in a tube and droplet spreading on a flat surface. The results show that the present approach leads to excellent convergence properties even with very coarse grids; furthermore, the results agree well with asymptotic analysis, with those obtained with a method for direct numerical simulations (wherein an adaptive grid is used) and also with experiments.