Modelling aerated flows with smoothed particle hydrodynamics

Modelling aerated flows is a complex application of computational fluid dynamics (CFD) since the interfaces between air and water change rapidly. In this work, the simulation of aerated flows with the smoothed particle hydrodynamics (SPH) method is investigated with a focus towards the application in engineering practice. To prove the accuracy of the method, the processes of air entrainment and rising air bubbles are studied. Through monitoring the evolution of the bubble contours it is shown that the novel approach of adding artificial repulsion forces at the interface does not alter the dynamics but stabilizes the flow. Building on these fundamental processes we extend the discussion to practical applications with a special focus on forced aeration. Since the employment of a detailed SPH model to practical problems remains out of bounds due to the high computational demand, we propose a combined experimental and numerical study where experimental bubble characteristics are imposed on the numerical simulation. Based on the data of the conducted bubble column experiment, the computational demand is significantly decreased such that the oxygen consumption due to biokinetic processes can be modelled. The future perspective is to apply SPH to urban water systems, e.g., for simulating detailed processes in wastewater treatment and sewer hydraulics. doi: 10.2166/hydro.2015.132 ://iwaponline.com/jh/article-pdf/17/4/493/388432/jh0170493.pdf Michael Meister (corresponding author) Wolfgang Rauch Unit of Environmental Engineering, University of Innsbruck, Technikerstrasse 13, 6020 Innsbruck, Austria E-mail: michael.meister@uibk.ac.at

based simulation modelling, smoothed particle hydrodynamics model, two-dimensional numerical simulation NOTATION Bo ¼ bond number (À) c ¼ speed of sound (cm s À1 ) D ¼ diameter of bubble column (cm) Δx o ¼ distance between orifices (cm) E i ¼ representative field variable at position of particle i (À) ε ¼ particle repulsion strengthness parameter (À) g ¼ gravitational acceleration (cm s À2 ) h ¼ smoothing length (cm) H ¼ water surface level (cm) m i ¼ particle mass (kg) P i ¼ particle pressure (Pa) ρ 0X ¼ reference density (kg m À3 ) ρ i ¼ particle density (kg m À3 ) q ¼ dimensionless kernel distance (À) Q ¼ aeration rate (cm 3 s À1 ) r i ¼ particle position vector in space (cm) R ¼ bubble radius (cm) Re ¼ Reynolds number (À) sd ¼ particle sampling distance (cm) t ¼ elapsed simulation time (s) T ¼ reference temperature ( W C) v i ¼ particle velocity (cm s À1 ) to additional constraints and relations which again make the numerical solution challenging (Chanson ). 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  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

PRINCIPLES
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 v, If subindices denote the value of a field variable at the position vector of particle i such that E i ¼ E(r i ), then with SPH this value is computed by weighting the influence of neighbouring particles by 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.

Governing equations
The governing equations can either be derived by integrating the fundamental interpolation formula by parts (see, e.g., Monaghan ) or from the Eckart Lagrangian by following a variational principle (Monaghan ). The latter derivation ensures that the SPH method reflects the symmetry properties of the Hamiltonian, e.g., linear and angular momentum are conserved. First, the differential form of the SPH continuity equation describes the rate of change of a reference particle's density by weighting the physical parameters of its neighbourhood Even though the subindex j includes all particles except for the reference particle, only neighbouring particles contribute due to the finite influence radius of the smoothing kernel. The double subindex abbreviates the difference between vectors such that v ij ¼ v i À v j except for the kernel where W ij ¼ W(r i À r j , h).
Second, particle acceleration is governed by the momentum equation which includes a multi-fluid stabilization term proposed by (Monaghan ) The first term approximates the pressure gradient and directly derives from the Lagrangian. The second term is attributed to the standard single fluid viscosity, although the averaging of density and viscosity is altered to support fluid phases with a large difference in viscosities. Conversely, the third term is not physically derivable, but empirically added to ensure stability of the method. Even though WCSPH is fundamentally multi-phase, for large density ratios such as for air and water the discontinuity in the density distribution causes an unstable phase interface which is initially repressed by adding an unphysically high surface tension (Colagrossi & Landrini ). Later corrections include the introduction of auxiliary functions, e.g., a particle number density (Hu & Adams ) or a particle volume distribution to avoid the discontinuity (Grenier et al. ). Whereas these algorithms require an extensive reformulation of the method, the present method is strongly related to standard SPH except for the R ij term in the momentum equation which stabilizes the phase interface by applying the original artificial repulsion force by Monaghan () between particles of different phases Based on a sensitivity analysis we, however, adapt the repulsion strengthness parameter to 7:75 Á 10 À2 since the higher value 8 Á 10 À2 proposed by Monaghan () 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.
Following the update of velocities, in a Lagrangian method like SPH, particles are advected with the local fluid velocity such that Third, the set of governing equations is closed by relating density to pressure by an equation of state of the form To limit the density fluctuations δρ=ρ to 1% as required to model an effectively incompressible fluid, the liquid's speed of sound is set to c X ¼ 10 max i∈X v i j j (Morris et al. 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 (c Y =c X ) ≈ 14 which results in small time steps. The particles' physical parameters for water as liquid and air as gaseous phase are summarized in Table 1.
Using the Wendland kernel of order 4 with the ratio Wall boundary conditions and viscous interaction bubbles are assigned a solid ellipsoidal shape and experimental rising velocities. In this study, the following shear viscous forces are substituted in the momentum equation To account for the correct viscous interaction it is important to ensure that the ellipsoidal bubbles carry sufficient layers of air particles such that the influence radius of the kernel is exceeded.

Implementation
The results of this work are produced with the in-house devel- 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 This test case is particularly suitable for demonstrating that the artificial repulsion force successfully stabilizes but otherwise does not affect the flow dynamics.

Breaking wave
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 Despite the physical importance of using a multi-fluid model, Table 2 illustrates the high computational demand associated with this simulation.

Rising air bubble
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 sd ¼ 5:8 Á 10 À3 cm in both phases, a circular bubble of radius R ¼ 2:3 Á 10 À1 cm is placed at the bottom of a water-filled tank with a width and height of 6R and 10R, respectively. The liquid's speed of sound is set to       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.

LIMITATIONS OF SPH FOR PRACTICAL APPLICATIONS
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.

Air inlet and treatment of free surfaces
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 neigh- 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  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 Q ¼ 2:7 cm 3 s À1 .
The bubble dynamics, which are captured by a high speed camera, are summarized in Table 4 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      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.

CONCLUSIONS AND FUTURE WORK
This work provides an overview of the SPH modelling of aer- cantly 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. ), where the high density ratio between air and water is avoided by separately calculating density and pressure for each phase.