A GPU-Compatible Scheme for Coupling SPH to ...

A GPU-Compatible Scheme for Coupling SPH to Process-Based Wastewater Treatment Models Michael Meister, Massoud Rezavand, Daniel Winkler, Wolfgang Rauch Unit of Environmental Engineering University of Innsbruck Innsbruck, Austria [email protected]

Abstract—In this work the weakly compressible SPH method is applied for modelling the hydrodynamics of a wastewater treatment (WWT) plant. In addition to the stirrer induced reactor mixing, a process-based WWT model is solved to account for the biological treatment of the sewage. Based on the principle of mapping pollutant concentrations to SPH particles, a GPUcompatible scheme for coupling SPH to process-based WWT models is developed. Since SPH particles are advected with the flow, these pollutant concentrations are transported along the treatment basin. The coupled model is validated for an intensely mixed treatment basin. After reducing the mixing intensity, the effects of the change in hydrodynamics on the biological treatment efficiency are analyzed.

I. INTRODUCTION In recent years the range of problems, which are modelled with SPH, have significantly increased. Despite that the research focus often remains on fundamental model development, which is usually validated against fundamental test cases like the example of a dam break (see e.g. [1]–[3]). As a typical SPH test case this problem involves highly dynamic flows, which are confined to a finite simulation domain. Focusing on problems like that, the importance of large spatial and temporal scales for simulating real-life engineering problems is easily underestimated. Whilst flows with large complex geometries are occasionally addressed in the literature (see e.g. [4]–[6]), there is a lack of simluations which involve time scales in the range of several minutes or hours. In this work the weakly compressible SPH method is applied for modelling the operation of wastewater treatment (WWT) plant. For this application the method’s simulation properties, e.g. the intrinsic support of free surfaces [7], multifluid flow and the exact treatment of advection [8], are beneficial. At the same time this is a prime example for a problem with large spatial and temporal scales. The sewage, which flows into the treatment plant, is subject to a characteristic diurnal variation. Hence, a minimum simulation time of 24 hours is required. In addition to the mechanical treatment of the sewage, the biological treatment of the contaminated wastewater should be considered. The underlying biokinetic processes also occur on time scales of hours. Therefore, it is of interest to simultaneously solve the hydrodynamics of the treatment basin and a process-based WWT model, which accounts for the biological treatment.

The principle idea of linking the SPH hydrodynamics to a process-based WWT model is reported in [9]. Pollutant concentrations of the process-based WWT model are mapped to SPH particles. Therefore, the so-called activated sludge model (ASM) [10] is solved for each SPH particle in addition to the hydrodynamics. The pollutant concentrations are transported along the treatment basin, because SPH particles are advected with the flow. In [11] this concept is extended by source and sink terms for the biokinetic model, which depend on the position of SPH particles within the reactor. Thereby, a detailed spatial distribution of the pollutant compounds throughout the treatment basin can be computed. In this work we address the coupling method between the SPH and ASM1 model in more detail and demonstrate the GPU compatibility of the scheme. A profound mathematical description of the method and a comprehensive introduction to biological WWT are included. In addition, extracts of [9] and [11] are adopted herein. II. METHODS In this paper the SPH method is coupled to the processbased WWT model ASM1. Firstly, we briefly introduce the employed formulation of weakly compressible SPH. Secondly, the process-based WWT model is described in detail. Thirdly, the coupling between the hydrodynamics and the biological WWT model is discussed. A. Weakly compressible SPH as hydrodynamic solver The applied method is weakly compressible SPH with a differential form of the continuity equation: X X dρi = ρi Vj (vi − vj )∇i Wij ≡ ρi Vj vij ∇i Wij (1) dt j j The discretized momentum equation includes pressure gradient and shear viscous forces contributions [12]: X p˜ij Fext (ri , vi , t) dvi =− ∇i Wij + dt m mi i j | {z } | {z } pressure gradient stirrer mixing (2) X vij + mi ν ρ˜ij ∂W/∂|rij | + g |rij | j | {z } |{z} shear viscosity

gravitation

11th international SPHERIC workshop

The abbreviations are p˜ij = (ρi pj +ρj pi )(Vi2 +Vj2 )(ρi +ρj )−1 and ρ˜ij = 2ρi ρj (ρi + ρj )−1 (ρ−2 + ρ−2 i j ). In addition to the gravitational acceleration, external forces resulting from the stirrer induced mixing are added to the momentum equation. Further details with regards to this force term and details about the applied corrective algorithms are discussed in section III. The ODE system described above is integrated in time with a second-order velocity-Verlet scheme. The update of velocities and positions by a ∆t/2 step is followed by a ∆t density integration. The timestep is concluded by a second ∆t/2 update of velocities and positions. The set of governing equations is closed by a direct relation between density and pressure, which is referred to as the equation of state: " # 7 ρi ρ0 c2i −1 (3) pi (ρi ) = 7 ρ0 The variables ρ0 = 1000 kg m−3 and c = 100 m s−1 denote the fluid’s rest density and the speed of sound. In this work we apply the concept of generalized wall boundary conditions [13] to model solid boundaries and use a Wendland kernel of order 4 with a smoothing length equal to the particle sampling distance. B. Activated sludge model as process-based wastewater treatment model The aim of this work is to apply the SPH method for a WWT simulation. State-of-the-art wastewater treatment requires multiple cleaning steps. The treatment cycle of the sewage begins with mechanical grate discharge and primary sedimentation. However, the sole separation of particulate material is not a sufficient treatment since soluble pollutant compounds are not removed with this step. Therefore, in modern WWT plants a biological treatment stage is added to achieve high purification rates. In addition to the hydrodynamics, process-based WWT models have to be solved for modelling biological WWT. Prior to discussing the coupling with SPH, the fundamentals of process-based wastewater treatment models are briefly reviewed. A more extensive introduction can be found in [9], [11] and [14]. In biological WWT a high purification rate of the sewage is obtained by bacteria which clean the contaminated water. Different species of bacteria utilize substances in the wastewater for their growth and maintenance. Since these bacteria reside in a sludge phase this natural process, which is adopted for technically treating the sewage, is referred to as activated sludge process. The key feature is to maintain the bacteria for a sufficiently long time in the technical system. Important pollutants in the sewage are ammonia or ammonium, which are degraded through a two-step process into molecular nitrogen, which degasses into the atmosphere at the surface. In this degradation process the intermittent pollutant compound nitrate is produced. Among others the requirements for an efficient operation of the biological treatment stage are a high level of hydrodynamic mixing in the treatment basin.

Munich, Germany, June, 14-16 2016

Depending on the type of the basin the mixing is either achieved by injecting pressurized air or driven by mechanical stirrers. In this work we focus on basins of the latter type. Sophisticated process-based WWT models account for several biokinetic processes, which describe the evolution of pollutant compounds and different species of bacteria. Despite that the mathematical structure of such models is simple. For each process an ODE with a specified reaction rate, e.g. a population growth or decay, describes the evolution of the biokinetic compounds. Since several processes influence each other, a system of ODEs has to be solved. The biological WWT models which are widely used are the activated sludge models (ASM), which were developed in series [10]. In this work the ASM1 model is implemented as the biokinetic model. This model defines 13 compounds Cm , each of which is controlled by the reaction rate Rm . Mathematically, reaction rates determine the rate of change in concentration over time and are derived from 8 underlying biokinetic processes. In the original ASM1 model the effects of the local hydrodynamic solution on the evolution of the biokinetic processes is neglected. The widely used procedure is to define each treatment basin to be fully mixed at all times. In large WWT plants a number of treatment reactors are connected in series and the flow rates in-between are specified. For the l−th reactor of volume V and influent rate Q the biokinetic compounds Cm are evolved according to the following equation [15]: Ql (t) dCl,m (t) = Rm + [Cin,l,m (t) − Cl,m (t)] dt Vl (t)

(4)

Thereby the subindices l and m ∈ {1, ..., 13} specify the treatment basin and the biokinetic compound. The variable Cin,l,m defines the m-th compound’s concentration in the influent to the l-th reactor. Using this concept detailed mixing effects are neglected and hence the spatial distribution of pollutant compounds in a treatment basin cannot be described [16]. By linking the detailed reactor hydrodynamics, which are modelled with SPH, to the ASM1 model, we conceptually improve the method and cure this deficiency. C. Coupling SPH and ASM1 With SPH as hydrodynamic model and ASM1 as processbased WWT model the question of how to couple these methods arises. Eventually, the effects of the hydrodynamics on the treatment efficiency should be characterized. Due to the aim of computing a detailed spatial distribution of the biokinetic concentrations throughout a treatment basin, a promising concept is to assign the compounds’ concentrations, which are defined in the ASM1 model, to SPH particles (see Fig. 1). Since the SPH particles are advected with the flow, the compounds’ concentrations are transported along the reactor. In particular, the influent of sewage with a high pollutant concentration and the effluent of treated wastewater can be characterized by SPH particles entering and exiting the domain. In this section the basic workflow for linking the SPH and ASM models is described. Furthermore, the


(municipal) wastewater influent internal influent from circulation pi =60.000 Pa

treatment basin


further treatment in next basin

VR1 (t) ≈ 150 m3 stirrer induced reactor mixing

pi = 0 Pa Spatially resolved evaluation of pollutant concentrations for each auxiliary cell

Nitrate pollutant

Soluble substrate concentration

Fig. 1. Overview of the treatment basin and the local computation of pollutant concentrations. A rectangular treatment basin with a dynamic municipal wastewater influent is modelled with SPH. Each SPH particle carries pollutant concentrations, which are dynamically evolved according to the ASM model. In the post-processing step these concentrations are evaluated at the position of each auxiliary evaluation cell.

benefits of mapping pollutant concentrations to SPH particles are explained within the context of GPU computing. Fig. 2 visualizes the workflow of a WWT simulation. The first step is the simulation of the reactor’s hydrodynamics, which requires the scene definition and variation of the wastewater influent as input. A standardized format for particle data is used to store the hydraulic states relevant for the ASM computation. As the effects of the biokinetic processes on the SPH model (e.g. accumulation of sludge) are currently neglected, the SPH simulation can be independently solved from the ASM simluation. The requirements for the SPH solver are the implementation of input and output interfaces so that CPU as well as GPU solvers can be employed. The second step is to use the particle simulation data as input for the ASM computation. Biokinetic evolutions of the pollutants are computed by providing the compounds’ concentrations of the inflow and reactor fluids. The result of the computation is stored in the same data format for each SPH particle. The ASM solver implementation allows to reuse or import existing SPH computations for testing various setups for the process-based WWT model. As previously explained, the pollutant concentrations of the ASM1 model are mapped to SPH particles such that the process-based biokinetic WWT model is solved for each particle. Even though this coupling method ensures that the concentrations are transported along the treatment basin, source and sink terms for the in- and effluent of the sewage have to be included in the model description. As a basic principle (4) is used. This equation describes the rate of change of concentration at the level of a fully mixed treatment basin. Source and sink terms for the hydraulics are added to the reaction rate of the ASM1 model, which is adopted unaltered. We now extend this concept to the level of individual SPH particles. The first modification is straightforward: The compound’s reactor concentration Cl,m in (4) is substituted by

a particle’s concentration Cm,i . The second modification concerns the hydraulic term. We define an inflow and outflow area at the top left and top right of the treatment basin [11]. These zones correspond to small fully mixed sections of the basin, where the sewage enters and the treated wastewater leaves the reactor. Correspondingly, whilst SPH particles reside within these zones, source and sink terms have to be considered for the evolution of the pollutant compounds. The extension of (4), which now describes the evolution of the concentration for an individual SPH particle denoted by the index i, is as follows: dCm,i (t) = Rm (t)+ dt  Q (t)C (t) in in,m  ,   Vinflow area (t) −Qin (t)Cm (t) + Voutflow area (t) ,    0,

if particle i in inflow area

(5)

if particle i in outflow area else (if particle i in interior)

According to (5) it is distinguished between 3 different cases. If the SPH particle is neither within the inflow nor outflow area of the basin, the concentration changes solely according to the reaction rate. In case that the particle is located in the in- or outflow area, source terms for the wastewater in- and effluent are considered in the process-based WWT model. In line with this procedure SPH particles flowing into a basin are initialized with the average concentrations of all particles within the inflow area. For each SPH particle (5) is integrated in time using the first-order forward Euler method with a constant timestep of 0.5 seconds. Since each SPH particle carries different concentrations, a spatial distribution of the compounds is computed. Thereby it is important to emphasize that the reaction rates are adopted unchanged from the ASM1 model. A self-evident future enhancement is to derive spatially dependent reaction rates Rm (ri , t), e.g. by extending the SPH solver to multifluid flows. Of particular interest for adopting local reaction rates for ASM1 is the modelling of pressurized aeration and the settling of sludge [17]. Following the computation of the hydrodynamics and the coupled process-based WWT model, the post-processing of the results should be discussed. Fig. 2 defines the workflow of the simulation and illustrates that the SPHASE GUI can

Fig. 2. Workflow of a WWT simulation using the SPHASE software, which includes an SPH and ASM model solver. The results of the hydrodynamics and the biokinetic processes are post-processed with the SPHASE GUI (Figure of [18]).


be employed as post-processing tool. In addition to analyzing the hydrodynamics, the evolution of the compounds’ concentrations over time can be plotted. The total concentration value of the treatment basin is computed by averaging over all SPH particles which are within the reactor. Even though each particle’s pollutant concentrations can be checked, it is impractical to analyze the ASM1 concentrations at the level of individual SPH particles. Therefore, the reactor is divided into ASM evaluation grid cells of constant size to analyze the spatial distribution of the compounds throughout the basin (see Fig. 1). It is important to note that these cells do not influence the hydraulic computation, but are only constructed for evaluating the spatial distribution. The biokinetic concentrations at the position of each auxiliary grid cell result from the smoothing over neighbouring SPH particles according to P ˜ j Cm,j (t)W (rk − rj , h) (6) Ck-th cell,m (t) = P ˜ W (rk − rj , h) j

Here rk denotes the centre point of the k-th grid cell. A ˜ which yields a cut-off radius of 6 neighsmoothing length h bours is specified. III. COMPUTATIONAL ASSETS AND CHALLENGES The simulation of a WWT plant is a good example for a problem with large spatial and temporal scales. Firstly, the underlying biokinetic WWT processes are slowly dynamic processes. The typical time scale where changes occur are in the range of several hours to days. Secondly, a municipal WWT plant is subject to a characteristic diurnal variation in influent contamination and flow rate, which is attributed to the humanly day-night rhythm. Therefore, a minimum of 24 hours of simulation time are required to account for biological WWT processes. Typical volumes of biological treatment basins are 100 m3 to approximately 800 m3 . A WWT plant consists of a number of such treatment basins which are connected in series. Consequently, the computational requirements for a WWT simulation require a trade-off in the level of detail of the SPH model. The compromise taken to tackle this challenge is to specify sampling distances in the range of 10 to 15 centimetres with around 15 neighbours per particle for a 2D simulation. This configuration corresponds to the lower resolution threshold for achieving representative reactor mixing for benchmarking the present coupling method. In favour of the coarse resolution we perform a 2D simulation and do not explicitly resolve the stirrer geometry with SPH particles, which would require a much higher resolution. Instead of that the effects of the stirrer are approximated by an external force field Fext (ri , vi , t) based on the principle of a PE controller [11]. In other words the mixing profile is already known from a detailed 3D mesh-based CFD simulation. A corresponding simplified flow profile is parameterized to dynamically control the external forces for each particle by comparing target to actual velocity.


A. Corrective algorithms for WCSPH A challenge for the application of SPH to WWT modelling is to maintain a regular particle distribution at the free surface over the entire simulation period of 24 hours. Even though the null pressure condition at the free surface is intrinsically satisfied [19], particles in proximity of the free surface have a lower number of neighbours which causes higher interpolation errors. Whilst the effects of these errors are negligible for highly dynamic processes occurring over short times, in a WWT simulation this dislodgement of particles can eventually lead to programme abortion. The regular mixing with bulk velocities in the range of 1 m s−1 is not sufficient to suppress this deficiency. In conjunction with increasing the speed of sound the application of a first order moving least squares (MLS) density renormalization [20] every 10 timesteps ensures that a regular particle distribution at the free surface is maintained over 24 hours. Even though higher order correction schemes further improve the accuracy of the free surface profile, the first order MLS well balances computational accuracy with efficiency. B. GPU implementation to increase computational efficiency Motivated by the high floating point arithmetic performance of modern GPU devices general-purpose computing on graphics processing units (GPGPU) is employed to solve highly parallel problems. This paradigm allows to harness the computing power of several trillion floating point operations per second (TFLOPS) for SPH simulations. The Nvidia CUDA platform provides a programming model and development environment that enables GPGPU development for Nvidia graphics cards and is used for the implementation of the SPHASE code. Although GPU accelerated SPH solvers exist, the discussed simplifications of the hydrodynamic model allow for specific optimizations of the implementation. CUDA GPUs implement a standard memory hierarchy that provides ample memory with relatively slow access times. For economic reasons the hierarchy is designed such that the speed of the memory increases the closer it is to the processor while the amount of storage decreases. The fastest memory levels are called caches and are usually controlled autonomously by dedicated hardware. The CUDA programming model makes an exception to this and allows the user to control a fast on-chip cache by software, called shared memory. For structured problems like matrix multiplications or mesh-based CFD methods this memory is often used to improve the performance of the GPU computation. In the computation of unstructured problems like SPH this cache is often neglected as the cache is either too small or the overhead of filling is too high. Due to the restriction to two dimensions and limited kernel support SPHASE is able to exploit this software controllable cache by employing a space filling curve for particle sorting. Fig. 3 shows a Hilbert space filling curve where every line is a connection of two points that have to be interpreted as grid cells containing particle data. Two particles have been chosen to visualize the concept of data locality, the center points


Fig. 3. Hilbert space filling curve. Every line is a connection of two points that have to be interpreted as grid cells, the solid blue line denotes the amount of data that fits into CUDA shared memory. Particles are visualized with bold points, the kernel support is indicated by the surrounding circle. (Figure of [24])

4 · 10−3

Time per iteration (s)

mark the particle abstraction into a grid cell, the surrounding circles visualize the relevant neighborhood. Neglecting the curve coloring and type, one can see that the two dimensional neighborhood is clustered into few coherent memory regions, which is beneficial for batch loading and caching performed by hardware caches. Interpreting the solid blue line in Fig. 3 as the amount of data that fits into software manageable cache visualizes the advantage of the approach. The light gray particle is calculated with all of the required data already cached and as such completely neglecting memory latencies. Whilst the black dashed circle does not exhibit this property, the ratio of particles that can be calculated from cached data may still be high. Due to the small number of particles in the simulation, the Verlet list (VL) concept [21], which is an a priori generated list of neighbors for each particle, is used. The higher memory consumption of this method is no limitation for the simulations, because the amount of simulated particles is very low to achieve a physical time of 24 hours. In [22] the performance of various algorithms for VL and cell linked list (CLL) is evaluated in favour of CLL due to better memory efficiency and similar performance for a CPU implementation. As GPUs benefit from iterating lists linearly [23] and memory restrictions do not apply to the present problem, the VL fits better to the architecture constraints than the nested loops in CLL implementations. A further optimization arises from the requirements of the application, namely that the integration of the biokinetic processes is performed every 0.5 seconds. For the example presented in section IV the time step resulting from the CFL condition is in the order of 3 × 10−4 and as such every ASM step requires approximately 1500 SPH iterations. Interim states of the hydrodynamic simulation are not required for the ASM1 computation and as such not transferred from the GPU device to the host memory. The dynamic time step size defined by the CFL condition prohibts to statically determine the exact simulation step. Transfer of the SPH data from the GPU to the host needs to be triggered by the host, which means that


SPHASE on GTX Titan X DualSPHysics on GTX Titan X SPHASE on GTX 750 DualSPHysics on GTX 750

2 · 10−3 1 · 10−3 5 · 10−4 2.5 · 10−4 1.25 · 10−4 104 Number of particles (-)

105

Fig. 4. Simulation performance chart. Performance comparison of SPHASE and DualSPHysics on an Nvidia GTX Titan X and GTX 750. The particle numbers result from the resolutions 0.02, 0.015, 0.01, 0.008, 0.006 and 0.005 m. (Figure of [24])

for each iteration the simulation time needs to be transferred and checked by the host. This results in a performance bottleneck for the implementation as the computation on the GPU requires less than 2 × 10−4 seconds. SPHASE estimates the approximate number of required iterations based on the CFL condition and the initial conditions for the sampling distance and the speed of sound. This approach allows to minimize the transfer of elapsed time while transferring a snapshot close to every 0.5 seconds.Standard optimization techniques for GPU programming have been implemented and are discussed in detail elsewhere [24]. Fig. 4 compares the average execution time for a single iteration at different resolutions for a dam break problem. The time per iteration is averaged for the simulation of 10 seconds of physical time. For low particle numbers SPHASE is more than 3 times faster than DualSPHysics using the GTX Titan X, although it is shown that this advantage constantly decreases. The lowest time per iteration for DualSPHysics on both devices is 5×10−4 , which means that SPHASE is more than 2 times faster when executed on the GTX 750 with scenes containing approximately 104 particles. It has to be noted that the results shown in Fig. 4 cannot be used to interpret the efficiency of the software implementation. While both solvers use the concept of WCSPH, DualSPHysics executes a simple Euler integration while SPHASE uses a second order velocity-Verlet scheme. The boundary conditions are conceptually similar but not identical. Furthermore, the 2D approximation of DualSPHysics causes a memory and computational overhead in comparison with SPHASE. The simulation of a WWT plant for 24 hours physical


time with a resolution of 173 × 39 particles is used to compare the performance of a parallelized CPU implementation with the discussed GPU implementation. The simulations are performed on a workstation equipped with 2 Intel(R) Xeon(R) CPU E5-2695 v2 @ 2.40GHz resulting in 24 cores (hyperthreading disabled), a GeForce GTX 750 and a GeForce GTX TITAN X graphics card. Using the CPU version of the code the simulation time for the low resolution equals 5.2 days, which rules out the computation of different scenarios. The GPU version requires 0.62 days for the same computation and shows superlinear scaling when increasing the resolution to 260 × 58 particles. While 2.24 times more particles are simulated, the wall clock time only increases by a factor of 1.35 to 0.84 days, which is due to unused resources of the GPU in case of the lower resolution. Extending the present simulations to 3D is unfeasible at the current state of research and technology. [4] use a high performance GPU (Nvidia GTX Titan) to simulate a tank of similar geometry, where 6 × 106 particles and 10 seconds of real time take 41 hours to compute. The resulting scale up factor is approximately 1.5 × 104 an thus even coarser resolutions will require years to compute simulations of 24 hours of physical time. IV. RESULTS: 24 HOUR WASTEWATER TREATMENT PLANT SIMULATION In this paper we present a WWT simulation of a rectangular treatment basin of 26 m length and 5.8 m height. A dynamic wastewater influent Qin (t) with a characteristic diurnal varia¯ in = 1800 m3 d−1 enters a tion and an average influent of Q treatment basin of volume V = 150 m3 . The process-based WWT model ASM1 and the hydrodynamics of the treatment plant are solved at two different resolutions: Firstly, 173 × 39 particles with a sampling distance of 15 cm are initialized on a square grid to characterize the treatment basin. This is the reference resolution, which unless otherwise stated, is used for generating the results. As depicted in Fig. 1 inflow particles enter the basin at the top left and exit at the top right of the reactor such that the overall particle number remains almost constant. Secondly, 260 × 58 particles with a sampling distance of 10 cm are specified. A. Model validation for strong mixing A brief reminder of biological WWT is that treatment is based on bacteria which utilize substances from the sewage. Therefore, a prerequisite for achieving a high treatment efficiency is to maintain an approximately uniform distribution of pollutant compounds and bacteria. In practice, among other methods, mechanical stirrers are employed for mixing the basins. As explained in section III the stirrer mixing is induced by an external force term which drives the flow. In [11] it is demonstrated that a single mechanical stirrer, which is placed in the centre of the treatment basin, thoroughly mixes the reactor and hence this configuration is adopted. For model validation the present method is compared with the classical ASM1 model’s results [10] which assume full


Fig. 5. Validation of the evolution of the compounds’ concentration. For a strongly mixed treatment basin (|vstirrer | = 1.5 m s−1 ) the pollutant concentrations in the treatment basin computed with the SPH based WWT model are compared to the classical ASM1 model, which assumes fully mixed reactor conditions (Figure of [11]).

mixing at all times. Using the SPH based WWT model this flow configuration can be approximated by specifying a stirrer bulk velocity of |vstirrer | = 1.5 m s−1 . Fig. 5 confirms that the predicted evolution of the biokinetic compounds well matches the reference results of the classical ASM1 model. The concentration of heterotrophic bacteria (XBH), which drive the biological treatment, is maintained and the soluble substrate (SS) as nutrient is sufficiently available. The compounds’ concentrations of ammonium (SNH) and nitrate (SNO) initially increase due to the high contamination in the influent, but these pollutant compounds are eventually degraded. B. Resolution independence of treatment efficiency Following the model validation the independence of the treatment efficiency on the particle resolution is checked. We repeat the WWT simulation using 260 × 58 particles and compare the pollutant compounds’ peak and final concentration after the 24 hour treatment cycle. Table I confirms that the predicted concentrations of ammonium and nitrate at t = 1 d are in good agreement. A relative deviation of 3.8 % is observed for the nitrate peak, which is within the limit of tolerance. Therefore it is confirmed that the reactor averaged overall treatment efficiency is independent of the two resolutions specified. At this point we underline that the higher particle numbers insignificantly increase the resolution of the hydrodynamics of the treatment basin, because the level of detail is mainly limited by the simplified description of the mechanical stirrer. In the next section it is demonstrated that a large number of SPH particles allows for the computation of a more detailed spatial distribution of the biokinetic compounds. C. Effects of reducing the mixing intensity So far the discussion focused on strong reactor mixing, which is driven by a stirrer with a bulk velocity |vstirrer | = 1.5 m s−1 . This high mixing intensity proved to be suitable for validating the present SPH based WWT model under



TABLE I R ESOLUTION INDEPENDENCE OF THE POLLUTANT CONCENTRATIONS . F OR A STRONGLY MIXED TREATMENT BASIN (|vstirrer | = 1.5 m s−1 ) THE OVERALL POLLUTANT CONCENTRATIONS IN THE TREATMENT BASIN DOES NOT DEPEND ON THE SPECIFIED SPH RESOLUTION .

Resolution 173 × 39 260 × 58 Relative deviation

Nitrate peak (g m−3 ) 5.53 5.75 3.8 %

)

S P H (v a ria b le m ix in g ) L in e a r fit

-3

n itra te p o llu ta n t p e a k c o n c e n tra tio n (g m

Nitrate at t = 1d (g m−3 ) 4.18 4.16 0.5 %

7 .7

7 .6

7 .5

7 .4

7 .3 0 .8

1 .0

s tirre r v e lo c ity

1 .2

1 .4

(m s -1)

Fig. 6. Dependence of the nitrate pollutant compound’s peak concentration on the mixing intensity. The maximum nitrate concentration at t = 0.71 d measured in the auxiliary evaluation cell positioned in the top centre of the treatment basin linearly decreases with higher stirrer mixing intensities (Figure of [11]).

the assumption of full hydrodynamic mixing. In practice, to reduce the costs of power consumption, the stirrer agitation should be reduced to the lowest level where a high treatment efficiency can be maintained. Using the present SPH based WWT model it is for the first time possible to analyze the effects of the mixing intensity on biological WWT processes. Herein we address the nitrate pollutant compound’s peak concentration in dependence of the mixing intensity as well as its spatial distribution throughout the treatment basin. A more comprehensive analysis of the results is found in [11]. Exemplarily, the peak concentration of the pollutant compound nitrate attained at t = 0.71 d in the top centre ASM evaluation cell of the treatment basin is analyzed. Fig. 6 confirms a linear decrease of the peak nitrate concentration with increased stirrer mixing intensities. The solid line approximates a linear function of slope −0.54 and with offset 8.08. The observed evolution is in accordance with the intuitive expectation that owing to the increased mixing intensity different concentration values within the reactor counterbalance such that the peak concentration decreases. The computation of the reactor averaged concentrations makes it possible to adjust the stirrer mixing velocities according to the target purification rate. The key innovation of the present method is that by reducing the mixing intensity to real-world values, a spatial distribution

Ammonium peak (g m−3 ) 12.6 12.57 0.2 %

Ammonium at t = 1d (g m−3 ) 6.06 6.05 0.2 %

of the concentration of biokinetic compounds arises. In the present work these concentration values are accurately computed for the first time, which is a significant improvement over state-of-the-art CFD modelling in WWT [25] and an exciting new application of SPH. Since each SPH particle carries concentration values, a detailed spatial distribution of the compounds is computed. Due to the large number of particles it is infeasible to analyze the treatment efficiency at the level of individual SPH particles. Using the SPHASE code as post-processing tool we visualize the pollutant concentrations throughout the treatment basin at the position of auxiliary grid cells (see Fig. 1). These concentrations are computed according to (6) by kernel-weighting over SPH particles which are located within each cell. Apart from quantifying peak concentrations it is of interest to locate local deficiencies in compounds’ concentrations. Since each SPH particle carries concentration values for the pollutant compounds, a detailed spatial distribution of the concentrations is computed. The method proposed for post-processing the results was introduced in section II-C. According to (6) the concentrations values are evaluated at the position of each auxiliary evaluation cell. Fig. 7 shows the spatially resolved concentration of nitrate after simulating a 24 hour diurnal variation with a stirrer mixing velocity of |vstirrer | = 0.7 m s−1 . The concentration values computed at the position of auxiliary cells form the basis for the contour plot shown. The results confirm that large parts of the basin are sufficiently mixed such that the pollutant is well degraded. The deviations at the edges of the basin are attributed to the reduced mixing in proximity to the walls. In these zones there is a reduced influent of wastewater with a high nitrate concentration and hence the degradation proceeds more quickly. Further investigation reveals that the magnitude of the nitrate concentration coincides with the reference values

4.2

2.5

Fig. 7. Spatial distribution of the nitrate pollutant compound. For a low mixing intensity of |vstirrer | = 0.7 m s−1 the spatial distribution of the pollutant concentration is evaluated at the position of the auxiliary evaluation cells in the basin (t = 1 d) (Figure of [11]).


for fully mixed basins. V. C ONCLUSION In this work SPH is coupled to the activated sludge model no. 1 (ASM1) for simulating a wastewater treatment (WWT) plant. The coupling method is based on mapping the pollutant concentrations, which are defined in the ASM1 model, to SPH particles. The ASM1 model is solved at the level of individual particles. Depending on a particle’s position within the reactor, source and sink terms are added to the model description. The coupling scheme is compatible with GPU computings and hence a 24 hour simulation time can be attained. The present WWT model is validated for an intensely mixed treatment basin, where detailed reference results are available. We then show that a reduced mixing intensity results in a linear increase of pollutant peak concentrations. In particular, we confirm that the pollutant concentrations can no longer be assumed to be fully mixed throughout the treatment basin. A detailed spatial concentration of the pollutant compounds is computed. Thereby, the treatment efficiency of a WWT plant can be closely monitored and problem areas within the basin can be identified. ACKNOWLEDGEMENT This research is part of projects which are funded by the Austrian Research Promotion Agency (FFG) [850738] and the Austrian Science Fund (FWF) [P26768-N28]. R EFERENCES [1] A. Crespo, M. Gómez-Gestaira, and R. Dalrymple. Modeling dam break behaviour over a wet bed by a SPH technique. Journal of Waterway Port Coastal and Ocean Engineering-ASCE, 134(6):313-320, 2008. [2] A. Ferrari, L. Fraccarollo, M. Dumbser, E. Toro, and A. Armanini. Three-dimensional flow evolution after a dam break. Journal of Fluid Mechanics, 663:456-477, 2010. [3] W. Jian, D. Liang, S. Shao, R. Chen, and X. Liu. SPH study of the evolution of water-water interfaces in dam break flows. Natural Hazards, 78:531-553, 2015. [4] A. J. Crespo, J. M. Dom´ınguez, B. D. Rogers, M. Gómez-Gesteira, S. Longshaw, R. Canelas, and O. Garc´ıa-Feal. DualSPHysics: Open-source parallel CFD solver based on Smoothed Particle Hydrodynamics (SPH). Computer Physics Communications, 187:204-216, 2015. [5] J. M. Dom´ınguez, A. J. Crespo, D. Valdez-Balderas, B. D. Rogers, and M. Gómez-Gesteira. New multi-GPU implementation for smoothed particle hydrodynamics on heterogeneous clusters. Computer Physics Communications, 184(8), 1848-1860, 2013. [6] A. Hérault, G. Bilotta, A. Vicari, E. Rustico, and C. Del Negro. Numerical simulation of lava flow using a GPU SPH model. Annals of Geophysics, 54(5):600-620, 2011. [7] J. Monaghan. Simulating free surface flows with SPH. Journal of Computational Physics, 110:399-406, 1994. [8] A. Colagrossi, and M. Landrini. Numerical simulation of interfacial flows by smoothed particle hydrodynamics. Journal of Computational Physics, 191(2):448-475, 2003. [9] M. Meister, and W. Rauch. Wastewater treatment modelling with smoothed particle hydrodynamics. Environmental Modelling and Software, 75:206-211, 2016. [10] M. Henze, W. Gujer, T. Mino, and M. van Loosdrecht (Eds.). Activated sludge models ASM1, ASM2, ASM2d and ASM3. Scientific and technical report No. 9., International Water Association, London, 2000. Publishing, London. [11] M. Meister, D. Winkler, and W. Rauch. A spatially resolved biokinetic wastewater treatment model based on smoothed particle hydrodynamics. Environmental Modelling and Software, submitted, 2016.


[12] X. Hu, and N. Adams. An incompressible multi-phase SPH method. Journal of Computational Physics, 227(1):264-278, 2007. [13] S. Adami, X. Hu, and N. Adams. A generalized wall boundary condition for smoothed particle hydrodynamics. Journal of Computational Physics, 231(21):7057-7075, 2012. [14] J. Makinia. Mathematical Modelling and Computer Simulation of Activated Sludge Systems. IWA Publishing, London, 2000. [15] O. Levenspiel. Chemical reaction engineering. Wiley, New York, 1972. [16] F. Hurtado, A. Kaiser, and B. Zamora. Fluid dynamic analysis of a continuous stirred tank reactor for technical optimization of wastewater digestion. Water Research, 71:282-293, 2015. [17] J. Habermacher, A. Benetti, N. Derlon, and E. Morgenroth. The effect of different aeration conditions in activated sludge - Side-stream system on sludge production, sludge degradation rates, active biomass and extracellular polymeric substance. Water Research, 85:46-56, 2015. [18] D. Winkler, M. Meister, and W. Rauch. SPHASE - Smoothed Particle Hydrodynamics in Wastewater Treatment. Proceedings of the World Environmental & Water Resources Congress, West Palm Beach, 2016. [19] J. Bonet, and T. Lok. Variational and momentum preservation aspects of SPH formulations. Computer Methods in Applied Mechanics and Engineering, 180:97-115, 1999. [20] M. Gomez-Gesteira, B., Rogers, R., Dalrymple, and A. Crespo. Stateof-the-art of classical SPH for free-surface flows. Journal of Hydraulic Research, 48:6-27, 2010. [21] L. Verlet. Computer experiments on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules. Physical review, 159(1):98-103, 1967. [22] J. M. Dom´ınguez, A. J. Crespo, M. Gómez-Gesteira, and J. C. Marongiu. Neighbour lists in smoothed particle hydrodynamics. International Journal for Numerical Methods in Fluids, 67(12):2026-2042, 2011. [23] NVIDIA corporation. Best practices guide, 2012. [24] D. Winkler, M. Meister, and W. Rauch. gpuSPHASE A shared memory caching implementation for 2D SPH using CUDA. Computer Physics Communications, in revision, 2016. [25] J. Laurent, R. Samstag, J. Ducoste, A. Griborio, I. Nopens, D. Batstone, J. Wicks, S. Saunders, and O. Potier. A protocol for the use of computational fluid dynamics as a supportive tool for wastewater treatment plant modelling. Water Science and Technology, 70(10):1575-1584, 2014.