Good modelling practice in applying computational fluid dynamics for WWTP modelling

Computational fluid dynamics (CFD) models have been extensively and successfully used for many decades in a wide range of engineering disciplines including, but not limited to, automotive engineering, chemical engineering, mechanical engineering and aerospace engineering. In recent years, there has been a steady increase in CFD modelling in the wastewater treatment (WWT) field, both in academia (universities and research institutes) and industry (consultants and manufacturers). This growth is due to increased availability of quality CFD models, increasing power and affordability of computing power, and increased importance of drivers for optimization of municipal and industrial WWT unit processes. CFD models are powerful tools that have become increasingly available with commercial packages providing graphical user interfaces to assist with model development, operation, and post processing. However, WWT unit processes involve a complex interaction of fluid, mechanical, biological, and chemical processes. Therefore, despite the improved prospects for application of CFD, these models still require significant expertise and experience to produce good quality and, hence, useful results.

A limitation noted by Nopens et al. (2012) is that there are currently few CFD modelling experts within the environmental educational and academic sectors. This lack of educational support can lead to software misuse by inexperienced users. Misuse leads to erroneous outcomes, which can result in bad designs or processes that result in suboptimal hydraulic design, contributing to distrust of CFD models in the WWT field. While there are several textbooks available on the general use of CFD modelling, currently there is no textbook or scientific report available that addresses practical applications of CFD specifically for WWT. Based on this lack of available guidelines, the International Water Association (IWA) has recently chartered a working group on CFD modelling to address issues of applying CFD to WWT engineering as well as training of CFD engineers. In addition, the relevance of precise details within WWT were historically not as significant when discharge limits were less stringent, energy input was less critical, and the facilities were sized for future conditions.

In contrast, other industries with a longer history of CFD model use and more stringent fluid flow requirements have established clear guidelines for quality modelling. Key references include Roache (1997), the 1998 American Institute of Aeronautics and Astronautics Guide to Verification and Validation of CFD Models (AIAA 1998), and the ERCOFTAC Best Practices Guidelines (Casey & Wintergerste 2000) for CFD modelling.

Several key choices need to be made at the outset of any CFD modelling problem to guide the direction of the modelling process. These choices include: (1) spatial dimensionality; two versus three dimensions; (2) time averaged versus time marching; (3) the number of phases and species; (4) material properties and their interactions when multiple species are present; and (5) if additional transport or transfer models are needed beyond what is already included in the CFD package.

In terms of spatial dimensionality, a three-dimensional (3D) problem may not substantially require more user effort to set up as opposed to a two-dimensional (2D) problem. However, it may have a significant impact on computational requirements due to the increased degrees of freedom. Hence, one should check whether the assumptions made when modelling in 2D are appropriate for the overall model objective. If this is the case, a 2D model might be able to provide the answer to the objective in a fraction of the time and cost required to develop and obtain the results in 3D. It should be noted that there are some opportunities to infer 3D outcomes from 2D geometry, for example with areas of axial symmetry in circular tanks or if a tank section is representative of the complete rectangular tank width. In these cases the inlet and outlet flow distribution should be uniform and 2D as well, as non-uniformity or rotating flow features will be lost in the solution with a 2D assumption.

Another consideration is time. For some applications, a steady-state solution of time averaged conditions is appropriate, particularly when boundary conditions are not changing. For systems with a water surface that is planar and steady, a rigid lid approximation may be used where the water surface is fixed at a known level, with no shear stress applied at the plane that would be an interface with air. However, if the flow is gradually or rapidly varied, a free surface model may be required. If boundary conditions or the solution is dynamic, a transient solution approach may be more appropriate. Although transient and free surface models are more computationally expensive and require appropriate initial conditions, they may be the only approach when solving a characteristic system which cannot be solved implicitly as a steady state.

A third important consideration is the number of species and phases required. Phases are states of matter including gases, liquids and solids, whereas species are different materials having a common state. WWT primary flows are typically liquid water that may include additional dissolved or fluid species, as well as solid and gas phases. Modelling of species can range from straightforward for passive, non-reacting materials, to very complex with reactions and leading to changes of phase state. Multiphase modelling is generally complex due to the phase interaction. Multispecies and multiphase modelling can be done with either a Lagrangian or Eulerian approach. A Lagrangian model individually tracks the species and/or phases as discrete elements within the bulk fluid, whereas an Eulerian approach treats the species and/or phases as a second (or more) fluid(s) in continuum (Crowe et al. 1998). The material properties for the species and phases being modelled will dictate the approach required.

Material properties can be chosen independently for all species and phases based on known values for steady thermal and non-reacting conditions, or set as an active function of other properties. Density and viscosity of the fluid is fundamentally dependent on factors such as solids concentration and temperature. Thermal effects can be included through the use of energy equations. The impact of solids on density is relatively straightforward to include in Newtonian flows using an additional transfer equation to couple the fluid and solid states in the solution. Recent research indicates that using a fixed density analysis of activated sludge tanks (i.e. assuming density of water) can dramatically over-predict velocities and the degree of mixing (Samstag et al. 2012). The impact of solids concentration on viscosity and density has been studied (De Clercq 2003; Ratkovich et al. 2013; Markis et al. 2014) and any analysis of sedimentation and mixing should account for these processes. The same may be said for the limit at which the liquid phase may be treated as a Newtonian fluid (Bridgeman 2012). For solids concentrations above 2%, the fluid becomes non-Newtonian where the physics are more complex and must be accounted for. Stress yield models have been generally applied in the literature (Wu & Chen 2007; Eshtiaghi et al. 2013; Ratkovich et al. 2013), but there is a lack of consensus, particularly at higher solids concentrations.

The Lagrangian particle tracking method is most commonly used in WWT modelling for grit and diffuse secondary species and phases that do not have a significant influence on the bulk flow. Lagrangian models are also commonly used for tracking passive organisms’ exposure to UV light intensity (Ho et al. 2011) when modelling UV disinfection. The Lagrangian approach can become cumbersome when modelling a significantly large number of particles is required (Crowe et al. 1998) or they have a significant influence on the bulk flow characteristics.

An Eulerian approach is generally used in WWT modelling for multispecies flows and for higher concentrations of secondary phases and/or when there is strong coupling between the phases. The first Eulerian approach is use of a transport model. The governing equations for fluid motion of the primary phase are solved, and the other phases and/or species follow the advective flow and are modelled as separate species or scalars with transport equations that predict local concentrations accounting for both advection and diffusion. Since most flows of importance in WWT are turbulent, diffusion is usually thought to be dominated by turbulence, rather than by molecular diffusion. When reactions are included, the transported species could be transformed to other phases of matter as solid flocs are formed, or gases are released. The interaction between momentum and mass conservation equations for the primary phase is complicated by making density and viscosity a function of scalar concentration, which could include temperature, salinity, or high biomass concentrations. The approach is computationally efficient as there are a limited number of additional equations to solve.

The second Eulerian approach is using a mixture model in which the phases are treated as a single continuous phase that is a blend of the discrete phases. The probability of the existence of each phase is based on the phase's volume fraction at a point in space. In this approach, a single differential equation is solved for continuity of the mixture. Also, a single differential equation is solved for the momentum of the mixture and tracks the motion of each phase relative to the centre of the mixture mass. This is handled by introducing the concept of a diffusion or drift velocity. A disadvantage of the standard mixture model approach for multiphase transport is that it relies upon input of characteristic sizes for calculation of phase coupling. This can be complicated in wastewater applications where a wide range of particle diameters are present and changing as they flow through a tank (Wicklein & Samstag 2009) due to mechanisms of flocculation and deflocculation. Furthermore, additional empirical functions have to be introduced to evaluate hindered and compression settling regimes that can occur, e.g. in secondary clarifiers.

A Eulerian-Eulerian approach is a third option where momentum and continuity equations are solved separately for each phase. The Eulerian-Eulerian approach has been thought to be preferable for liquid-gas flows (Amato & Wicks 2009) since it can incorporate a range of typical bubble properties, e.g. surface tension, agglomeration and break-up. However it does come at significant computational cost due to the number of equations that must be solved.

A final choice is whether additional models apart from the CFD are required. If conversion of suspended and dissolved species is occurring (and of interest to hydraulics, rheology, or the process itself), conversion terms need to be added to their respective scalar equations. If settling is to be modelled, a settling velocity in the direction of gravity should be introduced. If kinetics are to be modelled, source and sink terms are to be included in the respective scalar equations. These terms can easily be recovered from a stoichiometric matrix if the model is available in that form. It is the vertical summation of the products of stoichiometric terms and process rate expressions. Further details on kinetics can be found in the IWA STR No. 9 (Henze et al. 2000).

Once key problem elements have been identified and choices made, the CFD model can be built. CFD models are used to solve a numerical form of the fundamental equations of fluid flow known as the Navier–Stokes equations, which are based on momentum balances. For most practical engineering problems being addressed today, the Reynolds-averaged form of the Navier–Stokes (RANS) equations are used (De Clercq 2003; Khan et al. 2006; McCorquodale et al. 2007; Ho et al. 2011; Liu & Garcia 2011). These equations do not form a closed set (ASCE Task Committee 1988), owing to the non-linearity of the original Navier–Stokes equations and their temporal averaging. Current CFD models solve additional equations representing the turbulence characteristics of a flow field, which is a key parameter in determining the nature of flow, eddy formation, circulation, flow separation, and flow interaction with surfaces.

In their general form, the RANS equations cannot be solved analytically. For solution of the RANS equations, it is necessary to discretize the flow domain into a computational mesh to define the actual locations where equations of flow will be solved. CFD models approximate the differential equations by numerical methods, that resolve the equations into a set of algebraic equations or time-dependent ordinary differential equations (Lomax et al. 2003). These equations are solved to provide hydraulic information (e.g. velocity, water surface elevation, and pressure) at a finite number of discrete points within the flow domain. Most current CFD models use the finite volume method as this approach allows the use of an unstructured computational mesh; however finite difference and finite element methods have also been also used.

There is increased interest in more complete transient versions of Navier–Stokes equations to capture the dynamic nature of turbulence through use of detached eddy simulation (DES) and large eddy simulation (LES) as computational resources continue to improve and increasingly complex problems are being addressed. However, these more advanced models have not found their way into routine use to date and will not be discussed in any significant detail here (Guerts 2004; Spalart 2009).

A computational mesh is developed within the flow geometry. There are many available approaches to mesh generation and computational cell type. It is important that the mesh is sufficiently refined in the region of interest, and of sufficient quality overall to provide a converging solution. In areas of the domain where solution variables are likely to have steep gradients, for example, a finer mesh is recommended to reduce propagation of numerical errors and promote convergence. The task of mesh generation is generally accomplished through the use of mesh-generating software that allows for definition of the model geometry, computational cell size, and grid density, and provides tools for grid quality analysis. Unstructured computational grids are the most common type, as they allow the greatest flexibility in defining the model domain and meshing properties.

Although modern computers have ever-increasing capabilities, a balance is always needed between overall mesh resolution and computational efficiency to generate a solution in a reasonable amount of time. Mesh independence analysis is generally recommended to verify the solution and is not dependent upon the number of nodes used to solve a given problem, i.e. the model solution changes with reducing the computational cell size (and increasing the overall cell count). A coarse grid may have a lower level of inherent stability (and higher degree of numerical diffusion), while a fine grid will be numerically more stable, but will require increased computational effort, both because of the increased degrees of freedom, but also because of an increased model stiffness (different order of magnitude in model process rates). Increasing the fineness of the mesh also imposes the need for smaller time steps to maintain stability for transient simulations. Consideration must be given also to the minimum number of computational cells across flow regions to maintain accuracy of the numerical scheme. For example, five cells as a minimum are required for second-order accurate solutions. The mesh should be finer in areas of interest and where sharp gradients may occur (e.g. inlet, outlets, free surface, near agitators, etc.), and smoothly transition between finer and coarser mesh regions. For most wastewater applications the bulk flow region is of primary interest and this is where most wastewater process reactions typically occur. Therefore, the flow details within the boundary layer are usually not of primary interest and consequently should be sufficiently refined for use of wall functions to determine the transition from no-slip at a wall to the free stream flow.

It is recognized that modellers often do not run mesh independence studies for all problems due to resource constraints, particularly where previous or similar representative geometries have been analysed, or reasonable engineering judgment can be utilized to develop the model mesh. Thus, when reporting CFD results, the mesh type should always be noted as well as the overall cell count and a detailed image of the mesh should be provided. This is a prerequisite for judging the quality of the CFD exercise as mesh quality is critical for accurate solutions and subsequent interpretation of the results.

Once the mesh is developed, the model solver needs to be set up to calculate the CFD solution. In general, a double precision solver with second-order discretization should be used for all wastewater problems. Many of the species modelled occur in low concentrations, implying the need for high numerical accuracy. Details of pressure-velocity coupling scheme and equation discretion can be found in Patankar (1980), Versteeg & Malalasekera (1995), and Lomax et al. (2003).

WWT typically requires several phases: wastewater, suspended solids and gas phase in aerated zones. If biokinetic changes are to be modelled, dissolved species must also be introduced and could transition between phases. The model phases will need their physical properties applied, including density and viscosity as a minimum. Additional factors may also be required when supplemental energy (like mixing) is included in the solution.

Almost any CFD problem in wastewater flows will typically require a suitable approach to turbulence. A wide range of turbulence models are available and routinely used. The selected model should be appropriate for the problems in question. There are numerous references that discuss the details of the various models (Versteeg & Malalasekera 1995; Wilcox 1998; Ferziger & Peric 2002). In current practice, two equation models such as k–ɛ and k–ɛ are most commonly used for wastewater modelling, although both computationally simpler and more complex models are sometimes used in some cases. Table 1 displays different types of turbulence models with advantages and disadvantages in their performance. Launder & Spalding (1972) provide an early overview reference on turbulence models. The basic turbulence models commonly used in hydraulic engineering have been reviewed by the ASCE Task Committee (1988) and Rodi (2000).

The flow field computed by the CFD model is a direct function of the flow conditions applied at the domain boundaries, known as boundary conditions. Typical fluid flow boundaries include inlet boundaries, outlet boundaries, pressure boundaries, symmetry boundaries, and wall boundaries. Inlet boundaries provide a velocity or mass flow rate in the three vector components into or out of the model domain, as well as turbulence characteristics. Pressure boundaries have constant pressure and turbulence characteristics, and flow can move in or out of the domain. Outlet boundaries only allow flow to travel out of the domain, and may have no pressure constraints or turbulence constraints. Symmetry boundaries allow no vector component normal to the boundary. Wall boundaries are considered solid with no flow through the boundary. The wall boundary type can either be no slip (zero velocity) with a roughness component, or a slip wall with no roughness component. Typically, a wall function is used to approximate the transition from the zero velocity at the boundary through the boundary layer into the free stream. These wall functions model the effective drag from the roughness at the wall, without requiring the large number of computational elements required to resolve the flow field within the boundary layer.

For wall boundary conditions, the no-slip condition (i.e. velocities equal to zero at the wall) is applied to all solid surfaces. At very small distances near the solid wall, a viscous sub-layer exists and is followed by an intermediate layer and turbulent core. In the viscous sub-layer, the flow is influenced by viscous forces and does not depend on free stream turbulent parameters. The velocity in this sub-layer only depends on the distance normal to the wall, fluid density, viscosity, and the wall shear stress. Typically, the viscous sub-layer is too thin to discretize and is therefore not included as a significant region in turbulence models in WWT reactors. However, the intermediate region, which includes effects of both the viscous sub-layer and the turbulent core, is more significant in size and requires more care in predicting the velocity and turbulence.

The intermediate sub-layer is bridged by utilizing empirical wall functions to provide near-wall boundary conditions for the mean-flow and turbulence transport equations. The purpose of these empirical functions is to connect the wall shear stress to the dependent variables at the near-wall grid node. This grid node must lie outside this sub-layer and reside in the fully-turbulent zone. There are two types of wall function provided in commercial CFD codes: equilibrium log-law wall functions and non-equilibrium log-law wall functions. Non-equilibrium log-law wall functions should be used when there is turbulent transport of heat, and species at a reattachment point.

Reactions, energy, biokinetics and multiphase flows require additional boundary condition information. The additional parameters are typically added at the flow and wall boundary conditions and could include the wall and fluid temperature, species concentrations, or phase information.

The selection of whether to perform a steady-state or dynamic (transient, time-dependent) CFD simulation is a key stage in the process. A steady-state simulation uses an iterative scheme to progress to convergence, but persistent oscillations in the convergence residuals may indicate dynamic flow (transient) characteristic and the model may need to be performed as a time dependent simulation. Dynamic simulations also take much longer to compute, so the initial choice of a steady simulation is generally the most practical option, unless dynamic conditions are part of the problem objective, such as when the input variables are simply changing (e.g. diurnal flow, a cascading flow, or flows with changing salinity, temperature or influent composition). A dynamic simulation approach can also be used to converge to the equivalent steady-state solution over time. Applying transient simulations can be very useful for multiphase applications, particularly biokinetic modelling when the interaction between multiple phases is complex.

Numerical solution methods are typically solved iteratively, the convergence error being the difference between the calculated solution and exact solution (Wilcox 1998). Unfortunately, the model will typically not reach the exact solution; therefore metrics are needed to judge when the model is sufficiently converged. Commonly, the residual error, the sum of the differences between the currently iterated solution and the previously iterated solution for each cell and its nearest neighbours (or next-nearest neighbours in the case of second-order schemes) is the first measure. This difference should reduce over time, while satisfying the main criteria that the mass flow should remain continuous. Finally, the solution at discrete points should be stable (Wilcox 1998).

The core CFD model, if set up properly, should not require calibration to solve the fundamental fluid flow equations accurately for a single species single phase flow with appropriate boundary conditions and an adequate mesh, as the parameters used are based on physics or known properties, and hence are not contestable. However, there are uncertainties in many typical wastewater applications that may require some level of calibration. Submodels for scalar quantities (e.g. solids) often require empirical fits that are subject to many factors that cannot be isolated for a given problem without some local data. For example, settling velocities can be measured on site. Solids profiles can be measured under one condition and the model compared with the measured data. In some cases, modifications in the fit parameters for the empirical models can improve the prediction of local site performance of the model. In general, however, the focus of data-model comparison for CFD purposes is model validation rather than parameter estimation or calibration as it is preferable to measure parameters for empirical model development through dedicated controlled experiments.

Data collection for validation is both time consuming and resource intensive as it typically requires advanced measurement techniques. Velocity measurements such as Acoustic Doppler Velocimetry or Laser Doppler Velocimetry (LDV) can be used. These require specific skills and expertise to conduct (Vanrolleghem et al. 2006) and have their limitations, e.g. the studied reactor should be relatively small and both water and reactor should be transparent (at least one wall for 2D observation) for LDV to be applicable. Additionally, spatial profiling of certain species within the reactor (e.g. solids concentration or dissolved oxygen) and reactive tracer methods have been used (Gresch et al. 2010; Liotta et al. 2014). Caution should be exercised when using a visual inspection of 2D–3D graphs of experimental data and model predictions as this approach is subjective. A quantitative approach that makes a numerical comparison of model error using an error evaluation technique is preferred.

Validation by solids and dye profiling (tracer test) is also a relatively inexpensive and well-documented procedure and can be used for validation of solids transport modelling for many WWT applications (Bender & Crosby 1984). Other non-reactive tracers can be used, such as lithium chloride, which is relatively expensive and needs special analysis equipment (e.g. atomic absorption) to measure low concentrations (mg/L) but yield good results. There are two main kinds of methods to simulate residence time distributions (RTD) with CFD: Lagrangian particle tracking (Thyn et al. 1998; Stropky et al. 2007) and solving the transport equation for a passive tracer with a scalar representing the tracer concentration (Talvy et al. 2007; Zhang et al. 2007). The second method is sensitive to numerical diffusion, so care is required when it is used. However, if the subsequent goal is to implement biokinetics, the best numerical RTD approach is to use the scalar method, because simulations that include biokinetics also utilize the transport equation (see further), and it is important to create similar conditions for both CFD-RTD and the subsequent integrated CFD biokinetics model (Le Moullec et al. 2008).

Caution should be used with tracer testing as similar RTD could be produced for multiple internal vessel geometries. Therefore, these tests should be used in conjunction with other methods when possible. It is important to note that the hydrodynamics depends not only on the shape and size of reactors and physical properties of the liquid, but also the energy engaged in the process and the flow rates (Potier et al. 2005). Hence, reactor hydrodynamics strongly vary over the course of the day, since liquid flow rate fluctuates significantly.

To date, CFD models have primarily been used for design analysis and troubleshooting, mainly due to the computational effort. However, more frequent use of CFD and the accumulation of knowledge should increasingly allow for optimization of design and operation, even in shorter time frames. This has already been applied to secondary sedimentation (e.g. Griborio & McCorquodale 2006; Xanthos et al. 2011) and mixing (Samstag et al. 2012), but will also be true for other processes. Optimization can both be in terms of process performance (such as improved conversion through ‘smart’ mixing and improved removal of a particular species) as well as in energy reduction (e.g. smarter location and operation of impellers, aerators and sensors serving as controller inputs (Rehman et al. 2014b)).

CFD can be used to develop biokinetic models of biological treatment reactors based on fully realized velocity fields in the reactor tanks rather than on simplified models such as the tanks-in-series (TIS) modelling approach typically used in the IWA Activated Sludge Models (Hangos & Cameron 2001) and other commercial models in wide use today. Care must be taken when incorporating biokinetic models as some have been derived with assumed mixing conditions. In some cases, the biokinetics have been ‘calibrated’ by forcing a TIS model to match available data through arbitrarily adjusting the model parameters within the degrees of freedom, neglecting the fact that TIS may miss key hydraulic details in the system (Alvarado et al. 2012). Generally, biological species (dissolved oxygen, ammonia, nitrates etc.) can be introduced as user defined scalars in the solver. Biological reactions then provide the production and/or consumption of these species, which are used as source/sink terms in the transport equations for each species (Rehman et al. 2014a). In multiphase systems, mass transfer between the phases can also be incorporated either using the user defined scalars or species transport within the solver. Efforts to combine CFD with ASM models are ongoing as well (e.g. Le Moullec et al. 2010b, 2011; Sobremisana et al. 2011)

Finally, next to the use of CFD as a stand-alone tool in system optimization, a second route will be to use CFD as a means to build up knowledge which can then serve to improve the simpler models we are currently using. The first examples have been reported in the literature where CFD models are used to build somewhat more detailed mixing models called ‘compartmental models’ (compared with the current practice of tanks-in-series) (Le Moullec et al. 2010a, 2011; Alvarado et al. 2012). In this way CFD models clearly have a complementary and supporting role to play in WWT modelling and by no means will replace the traditional ones in the near future. On the contrary, they are an important tool to extend our knowledge further and exploit the power of mathematical modelling for WWT systems as recently corroborated by Laurent et al. (2014).

The potential for CFD in WWT applications is great. However, to date this potential has not been fully exploited. One of the reasons for this, next to lack of training for practitioners, is lack of guidance with respect to GMP for water and wastewater applications. GMP should cover the complete process from objective definition to problem solution and visualization of results. In this paper we have tried to outline briefly the principal features of GMP for WWT. This should allow modellers new to the application of CFD in the wastewater field to conduct CFD studies. This will increase trust in these models and allow the exploitation of their full potential in the future.

This work has been completed as part of the work of the IWA Working Group on CFD.