Wastewater treatment modelling with smoothed particle hydrodynamics

In this paper a novel hydrodynamic wastewater treatment (WWT) model based on smoothed particle hydrodynamics (SPH) is presented. The hydraulics of the wastewater treatment plant is modelled in detail with SPH. The SPH solver is coupled to the activated sludge model such that the influence on biokinetic processes is described. The key innovation of the present WWT model is that both the biokinetics and the wastewater hydraulics are simultaneously solved for non-steady flows. After validating the present method against the software ASIM 5, the capabilities are demonstrated for a full-scale treatment plant simulation. We investigate the stirrer and aeration induced mixing within the reactor compartments as well as the resulting concentrations of the biokinetic compounds. Following the establishment of a local coupling between the hydraulics and the biokinetics, the biokinetic concentrations within a treatment plant can be spatially resolved with a high resolution. © 2015 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Software availability Name of the Software: SPHASE e Smoothed Particle Hydrodynamics Activated Sludge Engine Developers: Gregor Burger, Michael Meister, Daniel Winkler Contact Address: Michael Meister, Unit of Environmental Engineering, University of Innsbruck, Technikerstrasse 13, 6020 Innsbruck, Austria, Michael.Meister@uibk.ac.at Availability: on application to the author Year first available: 2015 Hardware requirements: PC with multi-core processor, optional: CUDA compatible GPU Operating systems: Linux, Mac OS X, Windows Programming language: Cþþ User interface: command line, graphical user interface Subscripts i: (constant) reference particle index FD, computational fluid dyus stirred-tank reactor; GPU, uation; SO, dissolved oxygen ics; SPHASE, SPH Activated WT, wastewater treatment. . Meister), Wolfgang.Rauch@ r Ltd. This is an open access article j 1⁄4 {1,...,N}: neighbouring particle index l 1⁄4 {1,2}: reactor index m 1⁄4 {1,...,13}: biokinetic compound index


Introduction
A key step in modern wastewater treatment (WWT) is to biologically treat contaminated water by a sequence of biokinetic processes. Depending on the pollutant concentration the hydraulics of the treatment plant, which includes the supply of oxygen and the retention of activated sludge, has to be controlled. Whilst the biokinetics are described in detail by the well-established activated sludge models (ASM) summarized by Henze et al. (2000), the reactor hydraulics is coarsely approximated by a series of fully mixed reactors (Levenspiel, 1972). In some instances the complexity of the hydraulic model is increased by including effects of axial dispersion. A breakthrough in WWT hydraulics, however, is achieved by the application of computational fluid dynamics (CFD) methods, which allow for the detailed modelling of processes like the settling and transport of solids (Larsen, 1977). Since effects of multiphase flow (Gresch et al., 2011) and biokinetic models (Le Moullec et al., 2010a,b;Wang et al., 2010;Sobremisana et al., 2011) can be incorporated, CFD is a promising tool for wastewater treatment simulations. Particularly if conventional models fail e.g. in case of dynamic inflows and significant variations in the contamination, the strength of CFD methods is evident (Laurent et al., 2014).
Existing research on CFD in WWT is restricted to mesh-based methods. To account for the biokinetics an advection-dispersion equation is solved on the existing grid (cf. Fig. 1a). Due to the restriction of the hydrodynamics to steady flows the computational demand and complexity of the model are low, but in return the application area for steady flows in WWT is limited. The application of the meshless Lagrangian SPH method gives rise to the following advantages: Firstly, SPH intrinsically supports free surface (Monaghan, 2005) and multiphase flows (Colagrossi and Landrini, 2003). Secondly, since particles move along with the fluid advection is treated exactly. Hence if SPH particles are assigned biokinetic concentrations, these compounds are transported with the flow which is a prerequisite for locally coupling the hydraulics to biokinetic models. Thirdly, since each particle's physical variables are accessible at all times, the computation of the hydraulics and the biokinetics can be carried out simultaneously for unsteady flows as outlined in Fig. 1b.
In this paper a novel WWT software, where the hydraulics are modelled with SPH and directly coupled to the ASM1 model by Henze et al. (1987Henze et al. ( , 2000, is proposed. The wastewater hydraulics are computed with a multi-fluid SPH algorithm, whose accuracy is demonstrated by an analysis of the mixing in WWT-plant basins. Following the description of the method we briefly summarize the biokinetic ASM1 model to specify the coupling to SPH. Following the validation of the model against the software ASIM 5, the capabilities of the present method are demonstrated for a full-scale treatment plant simulation.

Methods
SPH is a particle-based discretization method, which can be applied to CFD problems by approximating a continuous fluid through a discrete set of particles. The evolution of these particles, which are assigned physical parameters like density r, pressure p and velocity v, is controlled by weighting the influence of its neighbourhood. The weighting by the SPH smoothing kernel W(r,h), in this work a Wendland kernel of order 4 (Hongbin and Xin, 2005) is used, is a key feature of SPH that simplifies the equations of fluid dynamics to a system of discretized ordinary differential equations (ODEs). Thereby r ¼ [r x ,r y ] denominates the particle position vector and h determines the width of the kernel. Exemplarily, the SPH approximation of a field variable A i ≡ A(r i ) is given by the fundamental SPH interpolation formula: In equation (1) it is to note that the computation of the field variable's value is subject to a kernel dependent error due to the substitution of the d-function by the smoothing kernel. In addition, the transition from the integral to the sum results in a discretization error. Based on equation (1) the SPH discretization of the governing equations can be derived as e.g. demonstrated by Monaghan (2005).

Introduction to SPH
SPH is the method of choice to model the treatment plant's wastewater hydraulics. Even though the present weakly compressible SPH (WCSPH) solver is directly applicable to 3D, a 2D simulation is performed to reach the time scale of days arising in biokinetic processes. The SPH discretization of the governing equations yields a closed set of ODEs. Firstly, the rate of change of a fluid particle's density is controlled by the continuity equation: These densities are reinitialized every 10 steps based on the moving least squares technique (see e.g. Gomez-Gesteira et al., 2010) to enhance the stability of the method. Particle velocities are then computed by the momentum equation: The derivation of the discretization of the pressure gradient term is found in (Hu and Adams, 2007). For the description of the viscous term F n , including shear and artificial viscosities, the reader is referred to (Monaghan and Gingold, 1983). In WCSPH the governing equations are closed by a direct relation between pressure Fig. 1. The hydraulic computation of a breaking wave modelled with a mesh-based CFD algorithm (a) is compared to the meshless SPH method (b). Coupled to the biokinetic ASM model the SPH method supports unsteady flows such that the scouring of soluble biokinetic compounds can be modelled. and density which is provided by the following equation of state: The reference densities and viscosities of water and air particles are set to their physical values at the base temperature of 20 + C. The liquid's speed of sound, which controls the size of the dynamic timestep, is set to 10 times the bulk velocity such that c ¼ 30 m s À1 . Consequently, density fluctuations dr=r are limited to 1%, which is required to model an effectively incompressible fluid.
An important asset for the application in WWT is that the SPH method intrinsically supports multiphase flows (Colagrossi and Landrini, 2003). The future aim is to explicitly model both the air and sludge phase with SPH to locally couple the multi-fluid hydraulics to biokinetic WWT models. Presently, the settling of sludge is not resolved in the SPH simulation. However, as visualized in Fig. 2 air is injected in reactor 2 of the treatment plant. Apart from supplying oxygen required for the treatment process, the air phase is required for mixing the basin. A drawback of SPH in case of high density ratios as with air and water is that an instability at the phase interface arises (Colagrossi and Landrini, 2003). Even though the numerical instability can be suppressed by corrective algorithms (e.g. Hu and Adams, 2007;Grenier et al., 2013;Monaghan, 2013), a large number of particles per air bubble are required. The implications on the modelling of aerated flows in WWT have previously been discussed (Meister and Rauch, 2015). Since the arising computational cost is incompatible with the time scale of WWT applications, we developed a novel algorithm on the basis of Ihmsen et al. (2011) where continuity and momentum equation are independently computed for each phase. Since it is out of scope to describe the algorithmic details in this short paper, the readers are referred to the original work (Bergmeister, 2015). In the context of the WWT simulation it is important to understand that the instability at the phase interface is removed instead of suppressed. In consequence a reduced simulation with 1e10 particles per air bubble can be specified such that the time scale of the biokinetic processes are within reach. Due to the coarse resolution specified, presently effects of turbulence are not resolved. Previous work on turbulence modelling with SPH (for an overview see Violeau and Issa, 2007) is typically restricted to fundamental problems and single-fluid flows.

Biokinetic model
Different species of bacteria utilize substances in the wastewater for their growth and maintenance. The activated sludge process facilitates this natural process for technical treatment of sewage.
The key feature is to maintain the bacteria for a sufficiently long time in the technical system. An important pollutant in the sewage is ammonia or ammonium (SNH), which is degraded through two processes. Firstly, the nitrification process biologically oxidizes ammonium to nitrite with a subsequent oxidation to nitrate. Secondly, the denitrification process stepwise reduces the nitrate to molecular nitrogen. Besides a sufficient concentration of bacteria and protozoa, each process requires specific environmental conditions. For nitrification a high soluble oxygen concentration (SO ! 2 g m À3 ) is required, whereas the denitrification process needs anoxic conditions, i.e. no dissolved oxygen. Whilst these 2 processes are suited to explaining the degradation of the compound SNH, in real applications several processes influence each other. In spite of that the mathematical structure remains simple. For each process an ODE with a specified reaction rate, e.g. simply a first-order decay, describes the evolution of the biokinetic compounds. Typically, a system of ODEs has to be solved since several processes simultaneously occur. For further background reading the textbook by Makinia (2010) is recommended for a step-by-step introduction to the field.
In practice, the sequence of biokinetic processes is described by the well-established activated sludge models summarized by Henze et al. (2000). The ASM1 model (Henze et al., 1987) defines 8 biokinetic processes and 13 compounds, each of which's rate of change is controlled by the reaction rate. In this paper the process rates and kinetic parameters of the original publication are adopted. In addition, we implement the corrections identified by Hauduc et al. (2010). The compounds' initial and influent concentrations are specified according to Table 1. In the original ASM1 model the effects of the local hydraulics on the biokinetics are neglected, whereas the global hydrodynamics are approximated by a series of connected continuous stirred-tank reactors (CSTRs) (Levenspiel, 1972). Within a reactor with volume V and an influent rate Q the biokinetic compounds C m are defined to be fully mixed at all times: dC l;m ðtÞ=dt ¼ Q l ðtÞ=V l ðtÞ Â C in;l;m ðtÞ À C l;m ðtÞ Ã þ R m Thereby the subindices l2{1,2} and m2{1,…,13} specify reactors and biokinetic compounds. The variable C influent,l,m defines the m-th compound's concentration in the influent to the l-th reactor. Whilst the principle of the model, i.e. a series of CSTRs, is adopted herein, the reactor hydraulics and the transition in  For the modelling of the clarifier we apply a simple point settler (see e.g. Makinia, 2010) and set the sludge retention time to 12 days. This sludge retention time determines the hydraulic residence time of the sludge within the treatment plant and is an input parameter in the biokinetic ASM1 model. The detailed dynamics of sludge transport and sedimentation can be included in the SPH model, but in this case the spatial resolution has to be significantly increased to resolve the correct dynamics.

Coupling ASM1 to SPH
The SPH time-step is dynamically controlled by a Courant-Friedrichs-Lewy condition and varies in between 2,10 À4 s and 8,10 À4 s, whereas for the biokinetic model a constant time-step of 8.64,10 À1 s is used. Since the SPH solver is executed much more often, the hydraulic computation in between two steps of the biokinetic computation is stored and later passed on for the coupling. For validation purposes we average the reactor flow rates Q l¼{1,2} and volumes V l¼{1,2} and establish a global coupling between SPH and ASM1 by inserting the computed hydraulics, i.e. the reactor flow rates and volumes, into equation (5). This ensures that the influence of the hydraulics on the evolution of the biokinetics is accounted for.
SPH particles carry identical reaction rates as specified in the original ASM1 model. Their concentrations only differ based on the initial level of pollution in the influent and the progress on treatment depending on the particle's retention time in the plant. In the future SPH will be locally coupled to biokinetic models such that the reaction rates of SPH particles spatially differ. The aim is to characterize the influence of aeration and sludge settling on the evolution of the biokinetics. Consequently, new areas of application like the modelling of localized oxygen or bacteria deficits emerge.

Implementation
The algorithms are implemented in a self-developed software titled SPHASE (SPH Activated Sludge Engine) whose key features are firstly released in this work. The coding structure of SPHASE is highly parallel by using OpenMP and can be run on platforms like Linux, Mac OS X and Windows. A GPU accelerated version of the code is subject to ongoing development and will be released as it becomes available. The distinguishing feature of SPHASE is that SPH is coupled to an ASM modelling routine. Computed results are instantly visualized in a graphical user interface such that the hydrodynamics and biokinetic processes can be monitored. In addition, the user can interactively control physical processes like excess sludge removal or wastewater composition.

Full-scale treatment plant simulation
The wastewater hydraulics of the biological treatment stage is modelled by 2 continuous stirred-tank reactors (CSTRs) which are characterized by oxygen setpoints of 0 and 2 g m À3 . As inflow condition the diurnal variation in the contamination and inflow rate of Rieger et al. (2012), which is generated with the influent software of Langergraber et al. (2007), is used. As shown in Fig. 2 the mixing in the anoxic reactor, where the biokinetic process of denitrification is predominant, is controlled by mechanical stirrers which are modelled by external forces, whereas in the aerobic reactor oxygen is supplied by rising air bubbles (Bergmeister, 2015). Table 2 summarizes the treatment plant's specifications. An internal recirculation of 2Q in and a return flow rate of Q in is specified. For the SPH computation a particle sampling distance of sd ¼ 15 cm and a ratio of h=sd ¼ 1 is specified. The wastewater influent at the topleft corner of reactor 1 (see Fig. 2) is realized using the buffer zone technique (Federico et al., 2012). The height of the influent, which is kept constant throughout the simulation, is defined by 2 SPH particles. The assigned inflow velocities derive from the inflow rate Rieger et al. (2012) and are dynamically adapted due to the diurnal variation. No special algorithm is required for the effluent since an inclined ramp at the top-right corner of reactor 2 allows for a free discharge.

Model validation
Prior to discussing the results the present WWT treatment model is validated with the well-established software ASIM 5 (Gujer and Larsen, 1995). The identical biokinetic ASM1 model is implemented in the SPHASE code and in ASIM 5. With regards to the hydraulics the ASIM 5 computation is based on two CSTRs which are connected in series. The dynamic wastewater influent is approximated by an hourly adjusted step function, which is the highest resolution that can be specified. In addition, ASIM 5 restricts the reactor volumes to constant values. Contrary, the SPH solver of the present model resolves the hydraulics in detail. In consequence, the dynamic wastewater influent is imposed on every SPH timestep such that the actual load curve is considered. Furthermore, the variation in the reactors' water-levels are included in the simulation.
In Fig. 3 the concentration profile of the key pollutant compound ammonium is compared with the ASIM 5 solution. The evolution of the concentration is in good agreement, whereby implications of the model simplifications in ASIM 5 can be identified. In particular the concentration profile in reactor 2 has rough edges due to the coarse resolution of the influent. In the present model's results, however, the anticipated continuous evolution occurs. The marginal deviation of the concentration profile in reactor 1 is explained by the fact that contrary to ASIM 5 the present WWT model considers variations in the water-level. Since a similar level of agreement is found for other biokinetic compounds, these illustrations are not included in this short paper.

Simulation results
For selected biokinetic compounds Fig. 4 depicts the concentration over one 24 h cycle. The wall-clock time taken to simulate this one day cycle on a system with 24 Intel(R) Xeon(R) CPU E5-2695 v2 processor @ 2.4 GHz with a cache size of 30 MB is 6.4 days. The application of GPU computing (Dominguez et al., 2013;Vacondio et al., 2014) allows for real-time simulation of the problem by using the SPHASE code with a single Geforce GTX Titan.
Due to the high initial values the concentration of autotrophic and heterotrophic bacteria weakly declines during the first cycle before attaining a steady state such that a high treatment efficiency is maintained. Since the process of nitrification is predominant in CSTR 2, the nitrate concentration increases but remains within the target range. Through pre-denitrification nitrate is then degraded to molecular nitrogen in CSTR 1. The effects of the dynamic wastewater influent on the biokinetic model, which are computed with the SPH solver, are clearly visible. The effluent concentrations after secondary treatment, which are summarized in Table 3 for the soluble substrate, nitrate and ammonium, confirm the build-up of the treatment efficiency and the effluent of sufficiently clean water.

Conclusions and future work
The proposed coupling between SPH and the biokinetic ASM1   In effluent (g m À3 ) 1.92 9.35 0.41 0 model represents a breakthrough in hydrodynamic wastewater treatment (WWT) modelling. The key step of the present approach is to couple the hydraulics to the biokinetic conversion processes. The present WWT model is validated against the software ASIM 5 and found to be in good agreement. The capabilities of the method are demonstrated for a full-scale treatment plant simulation with a dynamic wastewater influent and pollutant concentration. The full potential of the model is revealed if the effects of the local hydraulics on the biokinetic concentrations are included. Numerical experiments on the accumulation of oxygen and sludge have been made, but an efficient local two-way coupling between SPH and ASM1 needs to be found in the future.