Modelling aerated flows with smoothed particle hydrodynamics

= bond number (−)

c

= speed of sound (⁠⁠)

D

= diameter of bubble column (⁠⁠)

= distance between orifices (⁠⁠)

= representative field variable at position of particle i (−)

ɛ

= particle repulsion strengthness parameter (−)

g

= gravitational acceleration

h

= smoothing length

H

= water surface level

m_i

= particle mass

P_i

= particle pressure

= reference density

= particle density

= dimensionless kernel distance (−)

= aeration rate

r_i

= particle position vector in space

= bubble radius

= Reynolds number (−)

= particle sampling distance

= elapsed simulation time (s)

= reference temperature (°C)

v_i

= particle velocity

= particle volume

= kinematic viscosity

= smoothing kernel

Subscripts

= reference particle index

= neighbouring particle index

= water phase subscript

= air phase subscript

If air is entrained in water, the flow dynamics are governed by multi-phase gas-liquid fluid. Since both phases are immiscible with different densities, the overall flow becomes compressible. Aerated flows are frequently encountered in hydraulics, most importantly in connection with weirs, chutes and spillways – thus contributing to energy dissipation. In environmental engineering, forced aeration of water by means of air bubble entrainment is essential to provide the oxygen needed for biokinetic processes (see, e.g., Henze et al. 1987). In the field of wastewater treatment particularly, aeration is decisive for both performance and cost of the treatment process (Bischof et al. 1996).

Chanson (2013) provides a broad introduction to the hydraulics of aerated flows, emphasizing its importance to the field of engineering. At a micro-scale the equations of fluid motion can be separately derived for both phases and then coupled with each other. It is crucial to supply an explicit model to account for the evolution of the interphases between air and water. However, due to the high computational effort, such direct numerical simulations can only be obtained in a small number of cases. As an alternative, averaged forms of the conservation equations can be used, resulting in two-phase Navier–Stokes equations. However, this simplified form of the governing equations gives rise to additional constraints and relations which again make the numerical solution challenging (Chanson 2013). In this paper, we introduce the smoothed particle hydrodynamics (SPH) method to study these problems. The advantages and drawbacks of the SPH method for modelling aerated flows of engineering applications are discussed and a novel hybrid solution is proposed.

SPH, which was introduced by Gingold & Monaghan (1977) and Lucy (1977), is a fully Lagrangian meshless computational fluid dynamics (CFD) method initially developed for astrophysical simulations. Since its foundation SPH has been successfully applied to fluid mechanics to simulate, e.g., free surface flows (Monaghan 1994; Gomez-Gesteira et al. 2010), transport phenomena (Tartakovsky et al. 2007) and multi-phase problems (Colagrossi & Landrini 2003; Hu & Adams 2007; Grenier et al. 2009). Owing to its Lagrangian nature, SPH has several advantages in comparison with Eulerian CFD methods and other numerical methods like the particle method of characteristics (see, e.g., Hwang et al. 2013): for example, advection is accounted for exactly and mass, momentum and energy are conserved (Monaghan 2005). With SPH the flow dynamics are modelled by dividing the fluid into a discrete set of particles which are evolved based on a weighted influence of their neighbourhood. By specifying corresponding densities and viscosities for particles of different phases, the method intrinsically supports multi-fluid simulations in the case of comparable phase densities. However, a drawback encountered by Colagrossi & Landrini (2003) for high density ratios (as given for air and water) is that the phase interface becomes unstable. In this work, we develop a novel multi-fluid SPH algorithm resembling standard weakly compressible SPH (Monaghan 2013). The key deviation is to extend the momentum equation by an artificial repulsion force which is applied between particles of different phases. While this procedure is known to cure the instability at the phase interface, it has not previously been applied to aerated flows.

Instead of focussing on a single detailed small-scale study, the purpose of the paper is to provide a practically oriented overview of the capabilities and limitations of the SPH method in modelling aerated flows. We begin by summarizing the present multi-fluid weakly compressible SPH (WCSPH) algorithm used. The accuracy of the method is then demonstrated by modelling the fundamental processes of air entrainment and rising air bubbles. The correct evolution of the bubble contours proves that the artificial repulsion force stabilizes the phase interface while otherwise not affecting the flow dynamics. Building on these fundamental examples the challenges and extensions required for advancing to practical applications are discussed. Referring to forced aeration (e.g., in wastewater treatment), the key issue of modelling continuous air inlets is addressed. We then raise the case study of a bubble column simulation as an example to identify the lack of computational effectiveness as a drawback of the method. Even though improvements like graphics processing unit (GPU) computing enhance the computing time, the detailed modelling of the aeration-induced mixing and the interfacial oxygen transfer remains out of bounds if conventional multi-phase SPH algorithms are employed. As a novel method of resolution we propose a combined experimental and numerical study to account for the homogeneous flow regime of a bubble column. Notwithstanding its limitations, this hybrid approach significantly reduces the computational cost such that the simulation of aerated flows in engineering applications becomes reasonable.

The principle of the SPH method is to simulate the hydrodynamics by dividing continuous matter into a discrete set of points which are referred to as particles. Each particle carries physical parameters like density ⁠, pressure P, velocity ⁠, volume V and constant mass m. Even though these interpolation points are initially distributed on a rectangular grid, SPH is a meshless Lagrangian method where the observer follows the movement of individual particles instead of retaining a fixed coordinate system. The time evolution of the particles' physical parameters is governed by the continuum equations of fluid dynamics which are numerically solved by interpolating from neighbouring particles with a smoothing kernel ⁠. Thereby is the position vector in space and h, which controls the width of the kernel, is referred to as the smoothing length. Required kernel properties are convergence to the -function in the limit of vanishing smoothing length, symmetry and partition of unity if a Shepard filter correction is used.

This fundamental interpolation formula allows the governing equations to be expressed as a set of discrete differential equations whose derivatives of variables are solely computable by kernel operations. This allows for an efficient numerical solution of the equations in a highly parallel framework.

Based on a sensitivity analysis we, however, adapt the repulsion strengthness parameter to since the higher value proposed by Monaghan (2013) restricts the evolution of the bubble contours. Minor deviations do not affect the stability of the method but, e.g., alter the evolution of the contours of a rising bubble. Note that this correction, which is controlled by the phase density difference and the parameter ⁠, slightly increases the particles' pressure in the neighbourhood band of the phase interface.

To limit the density fluctuations to 1% as required to model an effectively incompressible fluid, the liquid's speed of sound is set to (Morris et al. 1997). For and if the values of other parameters are also replaced, an equivalent equation of state holds for the gaseous phase. Note that for air-water flow, the speed of sound ratio is which results in small time steps. The particles' physical parameters for water as liquid and air as gaseous phase are summarized in Table 1.

The results of this work are produced with the in-house developed SPH code titled SPHASE which abbreviates Smoothed Particle Hydrodynamics Activated Sludge Engine (SPHASE). The code is subject to open source licensing by the end of the research project. The name SPHASE is attributed to the aim of establishing SPH as a numerical CFD engine in the field of wastewater treatment by coupling it to the description of biokinetic processes. Still, the code is very general and allows application to a variety of problems, including aerated flows. The key feature of SPHASE is the real-time visualization of the computation, which is performed in a highly parallel framework. Open multi-processing (OpenMP) is the employed application programing interface that provides multithreading of C + + code on various operating systems including Linux, Mac OS X and Windows. Constructs for workload balancing, thread creation and synchronization are provided by OpenMP on all parallelizable code parts. In SPHASE these parts correspond to the majority of the code which is supported by the highly parallel structure of the SPH method. The particle neighbourhood search, which is the limiting structure for performance, is optimized by performing the radix sorting algorithm provided by the thrust library on a uniform grid which leads to a cost of (Merrill & Grimshaw 2011) with n denoting the number of particles. To fully exploit the advantage of the instant visualization, real-time user interaction is an important feature of this code. Physical properties of particles are accessible at all times and inflow rates are adjustable during the simulation. For visualization the cross-platform application QT is employed. The code further complies to the H5Part format standard for exporting simulation data. Finally, SPHASE can be independently used as a post-processing tool for SPH data available in this format. Throughout this work we reference simulation runtimes for an Intel(R) Xeon(R) CPU E5-2695 v2 processor @ 2.4 GHz × 24 with a cache size of 30 MB and 32 GB of DDR3-1600 RAM. The operating system is Arch Linux 64-bit.

Entrained air bubbles and pockets are an important feature of aerated flows and occur if the shear forces at the surface exceed the resisting surface tension. In fluid dynamics, the modelling of entrained air is a key issue since it significantly affects the flow due to its intense interaction with the liquid. In nature, entrained air enhances the interfacial oxygen transfer which is required for the natural self-purification and mixing of a water body. Even though aerated flows are subject to extensive studies, detailed CFD simulations are recommended to advance future understanding (Chanson 2013). However, due to the complex highly dynamic processes involved the modelling of entrained air is challenging with Eulerian CFD methods, whereas the present SPH algorithm is a powerful method to capture these multi-fluid interactions. This is demonstrated by two examples: first, the modelling of entrained air is illustrated by the breaking of a wave. We then study the rising bubble problem to prove that the present method allows for an accurate simulation of the bubble contours and the coalescence with the free surface. This test case is particularly suitable for demonstrating that the artificial repulsion force successfully stabilizes but otherwise does not affect the flow dynamics.

The standard example to illustrate the power of SPH in simulating highly dynamic processes is the dam break. Since this test case is thoroughly discussed in the literature (e.g., Dumbser 2011; Kao & Chang 2012) but not addressed using the simple multi-fluid SPH algorithm (Monaghan 2013), the fundamental simulation characteristics are not repeated herein. It is, however, important to demonstrate that the artificial repulsion force cures the numerical instability observed at high density ratios as with air and water (Colagrossi & Landrini 2003) such that the multi-fluid characteristics, e.g., the entrainment of air and the fluid's interaction with the surrounding air, are resolved in detail.

Every 25 time steps a Shepard filter correction is applied to ensure zeroth-order consistency of the kernel interpolation (see, e.g., Gomez-Gesteira et al. 2010). The geometric configurations are chosen as by Colagrossi & Landrini (2003), i.e., the water block's ratio of length to height is 2, whereas the height and width of the simulation domain, which is filled up with air particles, is 3 and 5.37 times the height of the water block. In the multi-phase simulation, the resolution is varied in between 2.10⁴ and 2.10⁵ particles to check convergence. Figure 1 illustrates the pressure and velocity distribution of the water phase for the highest resolution with a sampling distance of in both phases. As depicted in Figure 2, at the time instant after the reflection of the wave at the wall the computed phase interface is in good agreement with the inviscid boundary element method solution by Colagrossi & Landrini (2003). Because a significant amount of air is entrained by the breaking wave, a multi-phase model is required to capture the physics; otherwise the entrained void region instantly collapses. In addition, we observe a more uniform distribution of fluid particles if the air phase is explicitly accounted for. The current implementation ensures that the two fluids remain distinctly separated. In contrast to other multi-phase SPH algorithms as demonstrated by Mokos et al. (2014), no unphysical voids occur in the air phase regardless of the chosen resolution, even if the ratio is increased. Despite the physical importance of using a multi-fluid model, Table 2 illustrates the high computational demand associated with this simulation.

As an example of forced aeration occurring in environmental engineering, we study the rise of an air bubble in a stagnant liquid tank. Specified by particles with a uniform sampling distance in both phases, a circular bubble of radius is placed at the bottom of a water-filled tank with a width and height of and ⁠, respectively. The liquid's speed of sound is set to ⁠. This configuration yields a Reynolds number such that a turbulence model is not required. The Bond number confirms that buoyancy forces predominate surface tension effects and thus a corresponding model is not included in the simulation.

Figure 3 depicts the particle distribution, pressure and velocity field of the rising bubble at ⁠. The bubble contour (black solid line) is compared to the SPH solution by Grenier et al. (2009). We analysed the sensitivity of the deformation of the bubble contours on the artificial repulsion strengthness parameter which is recommended to be set to (Meister & Rauch 2014). In consequence, the stiffness of the phase interface is relaxed and hence the bubble is deformed appropriately, i.e., it attains a horseshoe shape which eventually splits into three segments. The optimized algorithm's evolution of the bubble contours as a function of the dimensionless time is plotted in Figure 4. Owing to the symmetry of the problem, only the left half-plane is shown. The current method's results agree with other SPH (Colagrossi & Landrini 2003; Grenier et al. 2009) and level set method data (Sussman et al. 1994; Wang et al. 2009). Minor deviations like the underestimation of the necking at the lower end of the horseshoe shape are identified.

Even though the present method accounts for the correct bubble dynamics within the liquid, the description of the air particles' movement fails as the free surface is approached. This is attributed to the weighting of the particle neighbourhood which becomes void if neighbours are missing. This difficulty is addressed in the next section, but we mention at this stage that the coalescence with the free surface is well described if the surface is explicitly filled up with air particles. To demonstrate this the height of the simulation domain is reduced so that the rising bubble approaches the surface at ⁠. The bubble then begins to interact with adjacent air particles. Single dispersed water particles, which penetrate into the air surface (see Figure 5(a)), are attracted to the liquid by the artificial repulsion forces and thus the different phases remain distinctly separated. The coalescence of the bubble alters the particle alignment at the surface, but a regular distribution without any singular particles or voids is maintained.

In the preceding section, we demonstrated the accuracy of the method in modelling the fundamental processes of air entrainment and rising air bubbles. The application to practical problems, however, is difficult due to the limited description of air inlets and the high computational demand which are summarized in the next section.

To advance to practical problems the presented single bubble study has to be extended to account for the modelling of continuous air inlets. Since the evolution of a particle's physical parameters is governed by a weighted influence of its neighbourhood, air inlet particles interact in both ways with fluid particles. In consequence, it is not possible to directly insert or delete particles and except for special cases even open boundary conditions cannot be accurately imposed (Ferrand et al. 2014). We extensively studied this problem with a novel inlet algorithm based on the buffer zone technique. If inlet particles are assigned constant velocities and the effective hydrostatic pressure such that the pressure gradient drives the inflow, then stable inflow configurations are attained. It is, however, unfeasible to impose the experimental inflow pressure, because the simulation blows up at an early stage. Adversely the bubble formation and detachment is unphysical. Even though the size of the developing bubbles agrees with experimental data, the detachment is deferred. Large numerical errors associated with this scheme cause a deformation of the symmetric bubble shape and the breakup into dispersed particles.

A physically correct inlet scheme can be realized by driving the inflow by a large enclosed pressurized chamber (Das & Das 2009), which is reviewed by Meister & Rauch (2014) using the present method. Despite the eligibility of this algorithm for academic test cases, this method cannot be applied to practical problems since up to 34,000 gas particles are required for modelling the detachment of a single bubble. With the current resources available, this computational overhead is not manageable except for a few small-scale studies. In addition, this method is not applicable to the modelling of a continuous inflow of air as, e.g., required for simulation of aeration-induced mixing in environmental engineering.

If a continuous inlet of air is attained, it is necessary to treat the degassing of air at the free surface unless it is explicitly described by air particles. In the latter case, as demonstrated, the coalescence of the bubble with the surface is accounted for correctly. In turn, this procedure is computationally expensive and only applicable to a finite simulation domain which is filled up with air particles. If treated as a true free surface, air particles have to be explicitly deleted, either based on a number of neighbours or a surface detection criterion (Ihmsen et al. 2011).

Despite the recent increase in computing power and the establishment of GPU computing (see, e.g., Herault et al. 2010; Dominguez et al. 2013), the SPH modelling of bubbly air-water flow is currently limited to single bubble studies (Szewc et al. 2013) and the merging of two bubbles if a lower resolution is chosen (Grenier et al. 2013). In consequence, the direct application of conventional multi-phase SPH to the modelling of bubble columns remains out of bounds for all practical problems. To allow for simulations in the field of environmental engineering a novel multi-scale approach is proposed where the bubble characteristics are determined by a detailed multi-fluid simulation (see section ‘Rising air bubble’) and assigned to a separate low resolution model. Building on this approach, we study the homogeneous flow regime of a bubble column with a combination of experiments and the numerical SPH method. Despite the simplified flow regime we should emphasize the complexity of this example since effects of wall-bounded turbulence and the arrangement of aerators significantly affect the flow dynamics (for an introduction see, e.g., Kantarci et al. (2005)).

Benefits of this hybrid SPH model are an increase of both particle sampling distance and time steps such that the simulation of the collective movement of 35 bubbles becomes considerably cheaper than resolving a single bubble in detail. In contrast to the air concentration-based model by Chanson (1996) the present algorithm resolves the flow field, which is a prerequisite for coupling SPH to biochemical processes where the bubble-induced mixing of the biokinetic compounds is important. With the approach taken it is for the first time possible to reach the timescales of biokinetic processes which range from seconds to days. Further, the experimental continuous air inlet configuration can be attained which is an unsolved problem for conventional SPH models. The practicality of the approach is demonstrated for the reference experiment below.

A reference experiment was conducted to determine the bubble velocities and to locate the homogeneous flow regime. Since the findings are passed on to the numerical simulation, the experimental flow configuration should resemble a two-dimensional model. Hence a sparger pipe with five equally spaced orifices with a diameter of is fixed at the bottom centre cross section of a three-dimensional cylindrical plexiglass column (see Figure 6) with the geometric specifications given in Table 3. This configuration limits wall effects and confines the bubble dynamics to the layer coplanar to the aerator. In the experiment, the aeration is driven by a membrane pump and its rate controlled by a flow-meter. To attain a pseudo-steady bubble configuration without any unsteady structures, i.e., a homogeneous bubbly flow regime, the aeration rate is set to ⁠.

The bubble dynamics, which are captured by a high speed camera, are summarized in Table 4. In this flow regime, bubbles of approximately uniform size and velocity are equally spaced throughout the domain and a terminal velocity of is reached. Even though the bubble movement is approximately two-dimensional, the ellipsoidal formed bubbles rotate along their axes with a minimum and maximum width between 1.13 cm and 1.64 cm, respectively.

Based on the experimental bubble characteristics, a two-dimensional numerical SPH simulation is performed in the cross section coplanar to the aerator. Despite an otherwise identical geometrical setup as in the reference experiment, three-dimensional bubble characteristics like rotation are not included in the numerical model. Since the experimental bubble distributions and velocities are homogeneous throughout the domain, in the numerical simulation bubbles maintain a solid ellipsoidal shape with a height of and a width of 1.4 cm (see Figure 7) with an assigned velocity from the experiment.

In comparison with the direct computation of the bubble dynamics, this combined study allows for modelling the air-induced mixing in a bubble column with 35 single bubbles at a much lower computational cost than required for resolving a single bubble (see Table 5). First, the high density ratio between air and water is numerically challenging for the SPH method because of the discontinuity in density at the phase interface (Colagrossi & Landrini 2003). Regardless of the choice of the corrective algorithm (e.g., Colagrossi & Landrini 2003; Hu & Adams 2007; Grenier et al. 2009), all solutions require a large number of air particles to characterize a single bubble. Since adaptive resolution is not yet established for multi-phase flow, the high level of detail required for curing the numerical instability is defined on the whole simulation domain even if it is otherwise not required. In contrast, if a kernel influence radius of 2 particles per direction is specified, then with the current method only 66 instead of 4,900 air particles per bubble have to be used.

Figure 8 indicates that the chosen particle spacing is sufficient for resolving the velocity field of the bubble column and for ensuring a stable long-term simulation. Second, in weakly compressible SPH the speed of sound in air should be 14 times higher than for water (Colagrossi & Landrini 2003) and hence the size of the time steps decreases accordingly. On the other hand, no change of the speed of sound is required if bubbles are modelled as solid air particles. In spite of that, Table 5 shows that the restrictive speed of sound of the bubble column simulation is comparable to the single bubble study since in the former case the bubble velocity is higher.

This work provides an overview of the SPH modelling of aerated flows. Results are produced with the self-developed code titled SPHASE which abbreviates Smoothed Particle Hydrodynamics Activated Sludge Engine. The software development is geared towards the future application in wastewater treatment but generally stimulates the future establishment of SPH as numerical CFD engine in environmental engineering. Real-time visualization of the computation is an important feature of the code and the corresponding program specifications are provided in this paper.

In multi-phase SPH modelling an instability occurs at the air-water phase interface which is cured by applying a weak artificial repulsion force between particles of different phases. The accuracy of the present method is demonstrated by modelling the fundamental processes of air entrainment and rising air bubbles. We have studied the impact of entrained air on the flow dynamics for the example of a breaking wave. In the rising bubble problem where effects of surface tension are neglected since buoyancy forces are predominant, the bubble contours are in accordance with level set method and reference SPH data. Hence unphysical effects of the artificial repulsion force can be ruled out. The modelling of continuous air inlets and the high computational demand associated with the method are identified as key limitations for establishing SPH in this field. Consequently, the simulation of permanent aeration and the resulting interfacial oxygen transfer remain out of bounds, despite an efficient code implementation and the use of multi-GPUs. Hence, in this paper, we propose a combined experimental and numerical study to model the homogeneous flow regime of a bubble column. Since the experimental bubble characteristics are imposed on the numerical simulation, the computational demand is significantly reduced such that the timescales of oxygen consumption due to biokinetic processes can be attained. In the future, this novel approach is expected to allow for the general application of SPH to environmental engineering problems, e.g., for modelling wastewater hydraulics. For a universal application the method is extended by a bubble coalescence and splitting model. A possible approach is seen in an algorithm developed for the animation of air bubbles (Ihmsen et al. 2011), where the high density ratio between air and water is avoided by separately calculating density and pressure for each phase.

This research is part of a project that is funded by the Austrian Science Fund (FWF): [P26768-N28].

Name of the software: SPHASE – Smoothed Particle Hydrodynamics Activated Sludge Engine

Developers: Gregor Burger, Michael Meister, Unit of Environmental Engineering, University of Innsbruck

Contact address: Michael Meister, Unit of Environmental Engineering, University of Innsbruck, Technikerstrasse 13, 6020 Innsbruck, Austria [email protected]

Availability: on application to the author

Year first available: 2014

Hardware requirements: PC with multi-core processor

Software requirements: Linux, Mac OS X, Windows XP/7

Programing language: C ++

Program size: 245K

User interface: command line, graphical user interface