A CFD strategy to retrofit an anaerobic digester to improve mixing performance in wastewater treatment

The wastewater industry is expected to face unprecedented challenges over future decades, with worldwide demand for clean water increasing by 50% and food demand by 30% by 2030 (WWAP 2012). The need to mitigate and adapt to climate change means that the link between wastewater and energy must be addressed (Dapelo et al. 2019). The role of anaerobic digestion of sludge from wastewater treatment is relevant here as, whilst methane-rich biogas is produced, mixing is necessary for the digestion process and is responsible for 17–73% of digester energy consumption (Owen 1982). However, recent experimental data (Kress et al. 2018) show that it is possible to reduce input mixing power without compromising nutrient distribution and process performance.

Over recent years, computational fluid dynamics (CFD) has established itself as a powerful tool to investigate mixing in anaerobic digestion. However, limitations persist. For example, the volume of published research relating to gas mixing (Vesvikar & Al-Dahhan 2005; Zhang et al. 2008; Wu 2010; Dapelo et al. 2015; Dapelo & Bridgeman 2018a, 2018b) is limited when compared to work on other types of mixing (viz., recirculated and impeller-driven), despite gas mixing being widely used. Furthermore, there is no consensus on a method to assess mixing quantitatively. A classical method based on evaluation of the average shear rate (Tchobanoglous et al. 2010) has been considered as inappropriate due to the likely coexistence of local shear rates of different orders of magnitude inside the same digester (Bridgeman 2012; Dapelo et al. 2015). Instead, Dapelo et al. (2015) proposed splitting the computational domain into zones where local shear rate values are expected to exhibit only modest variations. Vesvikar & Al-Dahhan (2005) and Hurtado et al. (2015) used the dead volume criterion based on the relative magnitude of the velocity field. Terashima et al. (2009) defined the uniformity index from the distribution of a passive scalar tracer. Finally, Dapelo & Bridgeman (2018a) introduced the concept of relative occupancy, a qualitative criterion based on the comparison of the range of values of a passive scalar tracer throughout the domain.

There is no universal model for multiphase flow. The Euler–Lagrange modelling technique is a natural choice for flows involving a dispersed phase, and has the advantage of needing a relatively small quantity of empirical data to close the momentum equations (Andersson et al. 2012). However, single particles being simulated separately means that computational expense may become prohibitively high when the number of particles grows (i.e., >10⁶, depending on the hardware). The Euler–Euler approach can reproduce a wide range of flows, but needs a quantity of empirical information to close the momentum equations. For this reason, the Euler–Euler technique should be considered when no other more specific models are available (Andersson et al. 2012). There has been a resistance from the scientific community towards adopting the Euler–Lagrange approach because of the perceived high number of particles potentially making the computational expense too high (Liu et al. 2017). Whilst this can be true for generic agitated vessels, which are the topic of the cited articles, it is not the case for the specific problem of gas-mixed anaerobic digesters, where the number of bubbles inside the digester at any timestep was shown to be below 10⁴ (Dapelo & Bridgeman 2018b).

No systematic comparison work between Euler–Lagrangian and Euler–Euler models for gas mixing in anaerobic digestion has been reported in literature. Only Wang et al. (2017) compared the results of Dapelo et al. (2015) to their own Euler–Euler population balance model (PBM) through a limited number of plots. The scope of the comparison was to demonstrate the applicability of the authors' model to a non-Newtonian liquid phase within a wider study on activated sludge/water, rather than comparing two different modelling approaches, which explains the limited range of the comparison. The Euler–Euler PBM work was reported to be no more accurate than the Euler–Lagrange. A possible reason for the gap in numerical accuracy may lie in an issue, which was present in the laboratory-scale simulation work reported by Dapelo et al. (2015), but not in Wang et al. (2017): as described below, the cells along the central axis had to be constructed with a much larger size than away from the centre. However, no comparison of the two models' numerical expense was reported. Despite the advantage of a possibly improved accuracy, Wang et al. (2017) proposed a steady-state only model, which precluded transient analysis and investigations on time-dependent mixing.

Wei et al. (2019) applied an Euler–Euler model to a full-scale digester mixed through vertical-hanging gas lances, and suggested that the reactor presented in Dapelo & Bridgeman (2018a, 2018b) was ‘essentially a bubble column with a bottom-mounted nozzle’, and hence ‘different’. The computational domain comprised up to 76,000 cells. The numerical work was undertaken in parallel by a two-core Intel i5-2740 for up to 10 days. Wei et al. (2019) concluded that their numerical results were consistent with Dapelo & Bridgeman (2018a, 2018b).

A further problem of the computational work undertaken so far is the lack of full-scale validations due to the technical problems arising from taking measurements on opaque and biochemically hazardous sludge on operating industrial plants. A compromise solution has been to perform validation at laboratory-scale, and then apply the validated model at full-scale.

CFD modelling of full-scale, gas-mixed digesters has unequivocally demonstrated that there is consistent space for mixing design improvement (Dapelo & Bridgeman 2018a, 2018b). However, previous work reported in the literature has been mainly limited to demonstrating the bare feasibility of such improvement, giving only sparse and fragmented indications on how to practically perform mixing improvement. Consequently, such work cannot directly help CFD consultants or professionals in the task of producing recommendations for water companies – at least, not beyond the simple indication that such work is possible.

This limitation has been addressed in the work reported in this paper. First, the model development and validation of Dapelo et al. (2015) is reviewed. Then, it is used to build a thorough CFD strategy to assess the mixing quality of a full-scale, gas-mixed anaerobic digester and produce retrofitting recommendations for mixing improvement: in particular, disparate considerations coming from literature (Dapelo & Bridgeman 2018a, 2018b) are integrated together into a coherent, structured and generalized protocol. The dependence of the results on different choices of parameters is assessed (Dapelo & Bridgeman 2018b). Quantitative and qualitative criteria to assess mixing are evaluated (Dapelo & Bridgeman 2018a). This analysis shows that the outcome of the CFD strategy is unaffected by difficult-to-measure parameters, and that more classical criteria for mixing assessment and sludge modelling should be reconsidered. Finally, this work is generalized to produce a flowchart to improve mixing in a generic anaerobic digester for wastewater treatment as a practical guidance tool for consultancy work.

A multiphase model for gas mixing in anaerobic digestion was introduced in Dapelo et al. (2015) and later applied in Dapelo & Bridgeman (2018a, 2018b). Flocculation and sedimentation were ignored as these processes occur over days if not years, whilst the time span of the modelled mixing was 20 min. Consequently, sludge was modelled as a single liquid phase obeying the incompressible Navier–Stokes equations, with non-Newtonian rheology dependent on total solid content (TS).

Further experimental results (Achkari-Begdouri & Goodrich 1992) showed that sludge density differed by no more than 1% from the reference value of 1,000 kg m⁻³ for a wide range of TS values, and hence, the latter value was assumed.

Biogas bubbles were modelled as pointwise Lagrangian particles. It is generally accepted that the pointwise particle approximation is accurate provided that particle size is much smaller than cell dimension (Andersson et al. 2012). This condition was respected in the full-scale simulations presented in this work, but not in the laboratory-scale validation work, as bubble size there is comparable to cell dimension. However, Sungkorn et al. (2011, 2012) showed that this requirement can be relaxed if the number of particles in the system at any timestep is below 10⁴. This latter condition applied in the validation work.

Direct observation (Dapelo 2016; Sindall et al. 2017) showed that, in the laboratory-scale, quasi-spherical bubbles rose well-distanced and non-interacting. For this reason, two-way coupling was considered, and coalescence and breakup were ignored. The force acting on a single particle was modelled as the sum of added mass, buoyancy, drag and lift force as in Dapelo et al. (2015). A drag model specific for gas bubbles rising throughout a liquid (Dewsbury et al. 1999) was adopted. A spherical particle approximation was assumed for simplicity.

Turbulence was modelled as follows. As bubbles were modelled as pointwise, the details on the wake were unresolved. Reynolds stress tensor was found to be an effective way to include the effect of the unresolved wake in the simulation results (Dapelo et al. 2015), and in particular, the Launder–Gibson Reynolds stress model (Gibson & Launder 1978) was adopted.

Transient laboratory-scale simulations (Dapelo et al. 2015) were performed on a 20 cm diameter, 13 cm height cylinder (Figure S1 in the Supplementary Material). Bubbles were injected from the centre of the bottom surface at gas flow rates of 2.05, 5.30 and 8.63 ml s⁻¹. The bubble effective diameter was observed to be between 7.01 and 7.94 mm. Simulations were performed for a total simulated time of 60 s. The computational domain comprised up to 2.3 million cells. Mesh independence was achieved through standard Grid Convergence Index (GCI) tests (Celik et al. 2008) provided that every cell containing a Lagrangian particle remained larger than the particle's nominal volume, hence, the necessity of keeping the cells along the cylinder axis unrefined.

The computational work was undertaken in parallel on three dual-processor 8-core 64-bit 2.2 GHz Intel Sandy Bridge E5-2660 worker nodes with 32 GB of memory, for a total of 48 nodes, for up to 48 hours.

Transient full-scale simulations (Dapelo & Bridgeman 2018a, 2018b) were performed on a simulated industrial digester (Figures S2 and S3 in the Supplementary Material and Table 2). The computational domain comprised up to 400,000 cells and consisted of a wedge of radians with periodic boundary conditions at the wedge faces.

To the authors' knowledge, no experimental procedure is available to measure the bubble size inside a full-scale digester; accordingly, a matrix of simulations with different bubble sizes (viz., 2, 6 and 10 cm diameter), was performed. The impact of this arbitrary choice is evaluated in the ‘Results and discussion’ section. In all the cases, the bubble size was to be smaller than the size of the cells containing particles, and hence no adjustment was taken in meshing. GCI tests were once more used to verify mesh independence. Biogas was injected through a series of 12 nozzles placed at the bottom of the tank, along a circular manifold of 1.83 m radius.

Each simulation was run in parallel on 12-core Intel Xeon E5-2690 v3 Haswell sockets running at 2.6 GHz with 128 GB RAM, for a total of 36 cores, for up to 36 hours.

Below, a number of mixing criteria existing in the literature are reported. Either they were used within the bespoke CFD strategy, or they were shown to be unsuitable to assess mixing in a gas-mixed digester, see the following sections.

For the average shear rateTchobanoglous et al. (2010) recommended a value of 50–80 throughout the computational domain for optimum mixing. As such, this criterion can be defined as quantitative and absolute: quantitative because it provides a number (viz., average shear rate); absolute because the degree of mixing is assessed through evaluating such number against a pre-decided interval only.

Dapelo & Bridgeman (2018b) proposed a relative variation of this criterion, with the degree of mixing being assessed by comparing the value of average shear rate to the outcome of a benchmark simulation instead of evaluating it against a fixed interval. Given the outcome of a simulation from a specific configuration, mixing is considered better than the benchmark if the average shear rate is larger than the benchmark's; worse otherwise.

Shear rate range's specific volume was proposed by Dapelo & Bridgeman (2018b) as a qualitative, relative criterion to support the conclusions of relative average shear rate. Four shear rate intervals were defined as very low (0–0.01 ⁠), low (0.01–0.1 ⁠), high (0.1–1 ⁠) and very high (1–⁠). Of course, other subdivisions are possible. For each interval, the proportion of volume where the local shear rate falls within the interval (‘specific volume’), is evaluated. A qualitative assessment of mixing is given by comparing each interval's specific volume to the outcome of a benchmark simulation: a given geometry is better than the benchmark if specific volume of the very low shear rate is smaller, and specific volume of the very high shear rate is bigger, than the benchmark; worse otherwise.

Dead volume (Vesvikar & Al-Dahhan 2005) is a quantitative, absolute mixing criterion. The portion of the computational domain where the velocity magnitude falls below 5% of the maximum velocity at a given time is labelled as dead volume. Mixing is defined as good if the value of dead volume is smaller than a pre-decided threshold.

The relative occupancy qualitative, relative mixing criterion considers the relative population of different logarithmic intervals for the scalar tracer (Dapelo & Bridgeman 2018a). Following its definition, the scalar tracer distribution becomes uniform when the relative population of one interval (but not the one corresponding to the lowest values of !) becomes predominant. Comparison to a benchmark case is required to assess mixing quality.

The strategy comprised two steps. The first consisted of finding an ideal value of input mixing power with the current design. A series of simulations were run where inlet flow rates were selected as fractions q of the value reported in Table 2. For each run, mixing was assessed through relative average shear rate and shear rate specific volume criteria, both with the case q = 1 as benchmark. The value of q corresponding to the ideal value for input mixing power was identified.

The second step consisted of altering the design of the biogas injection system, using the ideal value of input mixing power found in the first step. Three modifications of the original design (‘orig’) were proposed: a new nozzle series was placed on a concentric circle of 5.49 m radius (‘new’); biogas injection was switched between original and new nozzle series every 1 min (‘1 min’); switching occurred every 5 min (‘5 min’). Dead volume, UI and relative occupancy criteria were used to assess the mixing quality of each of the modifications.

If required, direct observation of the flow patterns may be used to corroborate the mixing criteria analysis.

Validation was performed through comparison of simulated flow patterns (Dapelo et al. 2015) with experimental measurements performed through particle image velocimetry (PIV) (Dapelo et al. 2015) and positron emission particle tracking (PEPT) (Sindall et al. 2017). Qualitative comparison of velocity flow patterns to experimental data showed a good level of agreement (Figure S4 in the Supplementary Material). The central high-velocity column on the PIV experimental flow patterns is attributed to optical aberration in the PIV measurements and treated as an artefact (Dapelo et al. 2015).

Shear rate was averaged over concentric cylindrical layers (Figure 1). The reduction in accuracy near the centre can be attributed to above-mentioned optical aberration and reduced data sampling in PEPT (Dapelo et al. 2015; Sindall et al. 2017). Furthermore, the Euler–Lagrangian model's constraints on particle size mean that the mesh was not optimal for accuracy around the centre. Despite these issues, the agreement between simulated and experimental data is good.

The relative average shear rate mixing criterion was applied to the ‘orig’ full-scale scenario. The value of the average shear rate as a function of q is plotted for all the sludge rheologies in Figure 2.

No relevant dependence on bubble size can be observed, the only exception being the highest flow rate runs for 6 cm diameter bubbles in Herschel–Bulkley (H-B) rheology. It is concluded that mixing predictions are not affected by bubble size (Dapelo & Bridgeman 2018b). Hence, a considerable modelling simplification, consisting of limiting further analysis to the less computationally expensive configuration (viz., 10 cm bubble diameter), can be safely adopted. For increasing values of mixing input power, Figure 5 shows an increase of average shear rate up to q around 0.5. For ⁠, the increase slows or stops, thus indicating 0.5 as the ideal value for q.

Application of shear rate range's specific volume (Figure 3) confirms both the observation regarding independence from bubble size and ideal value for q: the specific volume curves tend to stabilize, with very low shear rate dropping consistently, for ⁠. It is concluded that mixing input power could be halved (viz., ⁠) without significant alterations of mixing efficiency.

As a confirmation of the conclusions above, the velocity flow patterns were analysed for the example case of 10 cm bubble diameter, 5.4% TS (Figure 4). It can be seen that a vortex is formed as sludge is pushed upwards by the rising bubbles, it migrates towards the external radially as it approaches the surface, and finally sinks to the bottom of the tank and ultimately towards the nozzle. A branch of this vortex sweeps the sloping base. This branch is not present (or is very weak) for ⁠, but is clearly present in both and ⁠. Thus, organic material is able to reach the bottom of the tank in both and ⁠, but not in ⁠. This means that (i) should be discarded; and (ii) and produce similar mixing effects, and hence, should be preferred as it requires less energy consumption.

Throughout the following part of this work, the value of was adopted and the bubble size was set to 10 cm diameter to reduce computational expense.

An example plot of dead volume and maximum velocity over time is reported in Figure 5. Large and uncontrolled dead volume variations over time were observed. This is attributed to its dependence on maximum velocity, which is expected to undergo large fluctuations due to the transient nature of the simulations and the large local variations of the velocity field around the bubbles. For this reason, dead volume was considered as an inappropriate criterion to assess gas mixing in a transient context (Dapelo & Bridgeman 2018b).

Two sets of initial conditions for the tracer concentration were defined as illustrated in Figure 6. was set to zero everywhere throughout the computational domain, except for a small number of cells, where it was set to 1 (‘seeds’). Two seedings were chosen: in ‘seed 1’, the non-zero cells were located near the digester's sludge inlet in order to mimic new sludge entering the system; in ‘seed 2’, non-zero cells were distributed regularly throughout the computational domain to reproduce the evolution of sludge already being present inside the digester at the zero timestep.

Figure 6 and Table 3 report the results of the UI tests for both the seedings. Both the seedings show a clear advantage of the switched strategies over the non-switched for all the rheology models. Also, an advantage of the ‘1 min’ strategy over the ‘5 min’ is clearly identifiable. Figure 6 also shows an advantage of the ‘new’ strategy over the ‘orig’; however, this is irrelevant for the purpose of mixing optimization because such advantage is completely overridden by the switched strategies' performance.

Figure 7 reports the outcome of the relative occupancy test. The observations offer qualitative support for the outcome of the uniformity tests. In addition, it was shown that Newtonian simulations display a clearly different behaviour in the evolution of the poorly mixed interval (black) from all the other (non-Newtonian) runs. Whilst it is well-known that non-Newtonian and Newtonian rheologies produce different flow patterns (Wu & Chen 2008), it was claimed in more recent research (Hurtado et al. 2015) that ‘there are not significant deviations between the results obtained between Newtonian and Newtonian models.’ However, the considerations above describe a clearly different behaviour in terms of mixing, and for this reason, the use of Newtonian models to reproduce sludge behaviour is discouraged.

The considerations above can be generalized into a flowchart, which can be applied to any gas-mixed digester and is readily usable by professionals and consultants (Figure 8). The mixing improvement strategies for the specific case presented within this work were limited to changing mixing geometry, but a general case may include alterations of a digester's geometry (e.g., baffles). The flowchart in Figure 8 referred to this general case as ‘altered geometry/injection’.

CFD modelling allows simple and practical solutions to reduce mixing input power considerably and increase process efficiency.

In order to generate reliable results, it is important to select carefully the mixing criteria and robustly reproduce the sludge's non-Newtonian behaviour.

Dead volume is an inappropriate method for determining the effectiveness of mixing. Instead the UI and relative occupancy concepts offer improved analysis.

A thorough CFD strategy to improve mixing performance was elaborated and described.

A general flowchart to improve mixing in a generic digester was provided. This flowchart can be readily used by professionals and consultants to produce recommendations for water utilities, operators and designers

The Supplementary Material for this paper is available online at https://dx.doi.org/10.2166/wst.2020.086.