A CFD strategy to retro ﬁ t an anaerobic digester to improve mixing performance in wastewater treatment

To date, mixing design practice in anaerobic digestion has focussed on biogas production, but no adequate consideration has been given to energy ef ﬁ ciency. A coherent, comprehensive and generalized strategy based on CFD modelling is proposed to improve mixing ef ﬁ ciency of a full-scale, uncon ﬁ ned gas-mixed digester for wastewater treatment. The model consists of an Euler-Lagrange (EL) model where biogas bubbles are modelled as the Eulerian dispersed phase, and non-Newtonian sludge as the Lagrangian continuous phase. Robustness tests show that mixing predictions are independent of bubble size. The CFD strategy comprises the assessment of different mixing geometries and a range of input gas ﬂ ow rates. Quantitative results show that simple retro ﬁ tting measures are able to achieve a signi ﬁ cant improvement in the degree of mixing with reduced mixing times, and consequently recommendations for best mixing geometry and gas ﬂ ow rate are given. A generalization to a generic digester is discussed in a form that is readily usable by professionists and consultants.


INTRODUCTION
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 ). The need to mitigate and adapt to climate change means that the link between wastewater and energy must be addressed (Dapelo et al. ). 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 ). However, recent experimental data (Kress et al. ) 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 ; Zhang et al. ; Wu ; Dapelo et al. ; Dapelo & Bridgeman a, b) 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. ) has been considered as inappropriate due to the likely coexistence of local shear rates of different orders of magnitude inside the same digester (Bridgeman ; Dapelo et al. ). Instead, (Dapelo et al. ) proposed splitting the computational domain into zones where local shear rate values are expected to exhibit only modest variations. (Vesvikar & Al-Dahhan ; Hurtado et al. ) used the dead volume criterion based on the relative magnitude of the velocity field. (Terashima et al. ) defined the uniformity index (UI) from the distribution of a passive scalar tracer. Finally, (Dapelo & Bridgeman a) 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. ). However, single particles being simulated separately means that computational expense may become prohibitively high when the number of particles grows (i.e., >10 6 , 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. ). 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. ). 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 4 (Dapelo & Bridgeman b).
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. ) compared the results of (Dapelo et al. ) 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 lab-scale simulation work reported by (Dapelo et al. ), but not in (Wang et al. ): 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. ) proposed a steady-state only model, which precluded transient analysis and investigations on time-dependent mixing. (Wei et al. ) applied an Euler-Euler model to a fullscale digester mixed through vertical-hanging gas lances, and suggested that the reactor presented in (Dapelo & Bridgeman a, b) 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. ) concluded that their numerical results were consistent with (Dapelo & Bridgeman a, b).
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 lab-scale, and then applying 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 a, b). 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 companiesat 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. ) is reviewed. Then, it is used to build a thorough CFD strategy to assess the mixing quality of a fullscale, gas-mixed anaerobic digester and produce retrofitting recommendations for mixing improvement: in particular, disparate considerations coming from literature (Dapelo & Bridgeman a, b) are integrated together into a coherent, structured and generalized protocol. The dependence of the results on different choices of parameters is assessed (Dapelo & Bridgeman b). Quantitative and qualitative criteria to assess mixing are evaluated (Dapelo & Bridgeman a). This shows analysis 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.

Multiphase model
A multiphase model for gas mixing in anaerobic digestion was introduced in (Dapelo et al. ) and later applied in (Dapelo & Bridgeman a, b). 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).
Sludge rheology was modelled through power-law: and Herschel-Bulkley models: where μ is the apparent viscosity, K the consistency coefficient, n the power-law index (n < 1 for sludge) and S the shear rate. In the Herschel-Bulkley model, flow occurs only if shear rate exceeds the critical value τ 0 . The values for the above-mentioned coefficients were defined following measurements presented in the literature (Landry et al. ; Baudez et al. ); furthermore, a Newtonian model was introduced as a comparison as some published work has chosen to ignore the non-Newtonian nature of sludge (Table 1). For the lab-scale validation work, a power-law transparent sludge substitute (carboxymethyl cellulose, CMC) was used in place of sludge to allow optical measurements. Different CMC concentrations were chosen to cover the range of the power-law parameter reported in literature (Table 1). Further experimental results (Achkari-Begdouri & Goodrich ) showed that sludge density differed by no more than 1% from the reference value of 1,000 kg/m 3 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. ). This condition was respected in the full-scale simulations presented in this work, but not in the lab-scale validation work, as bubble size there is comparable to cell dimension. However, (Sungkorn et al. , ) showed that this requirement can be relaxed if the number of particles in the system at any timestep is below 10 4 . This latter condition applied in the validation work.
Direct observation (Dapelo ; Sindall et al. ) showed that, in the lab-scale, quasi-spherical bubbles rose well-distanced and non-interacting. For this reason, twoway coupling was considered, and coalescence-breakup were ignored. The force acting on a single particle was modelled as the sum of added mass, buoyancy, drag and lift force as per in (Dapelo et al. ). A drag model specific for gas bubbles rising throughout a liquid (Dewsbury et al. ) 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. ), and in particular, the Launder-Gibson Reynolds stress model (Gibson & Launder ) was adopted.

Laboratory-scale meshing
Transient lab-scale simulations (Dapelo et al. ) 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. ) 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  (Dapelo & Bridgeman 2018b), and rheological properties of CMC solutions used for lab-scale work (Dapelo et al. 2015) TS/conc.  Table 2). The computational domain comprised up to 400,000 cells and consisted of a wedge of π=6 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' 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 taking 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.

Mixing criteria used within the CFD strategy
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 gas-mixed digester in the following Sections.
Average shear rate. (Tchobanoglous et al. ) recommended a value of average shear rate of 50-80 s À1 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 b) 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 b) 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 s À1 ), low (0.01-0.1 s À1 ), high (0.1-1 s À1 ) and very high (1-∞ s À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 very low shear rate's specific volume is smaller, and very high's specific volume bigger, than the benchmark; worse otherwise.
Dead volume (Vesvikar & Al-Dahhan ) is a quantitative, absolute mixing criterion. The portion of the computational domain where the velocity magnitude falls below 5% of the maximum velocity as 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 Uniformity Index (UI) is a quantitative, absolute criterion introduced by (Terashima et al. ). A scalar, numerical tracer χ is defined as obeying a non-diffusive advection-diffusion equation: whereũ is the velocity field. The uniformity index is then defined as: where χ i and V i are respectively the value of χ in, and the volume of, the i-th cell, χ is the average value of χ and throughout the domain and V the volume of the computational domain. As a consequence of the definition in Equation (4), the Uniformity Index spans between the values of 0 for total mixing and 1 for total inhomogeneity   Mixing is defined as good if the value of the uniformity index falls below a pre-decided threshold: accordingly, (Terashima et al. ) defined the thresholds of 90% mixing (UI < 0:10), and 95% mixing (UI < 0:05). The relative occupancy qualitative, relative mixing criterion considers the relative population of different logarithmic intervals for the scalar tracer χ (Dapelo & Bridgeman a). 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.

Mixing-improvement CFD strategy
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 was run where inlet flow rates were selected as fractions q of the Q max 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, uniformity index 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. 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. ; Sindall et al. ). 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.

Laboratory-scale validation
Optimization strategy first step: finding the ideal value of input mixing power 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 rheology. It is concluded that mixing predictions are not affected by bubble size (Dapelo & Bridgeman b). Hence, a considerable modelling simplification, consisting of limiting further analysis to the less computationallyexpensive 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 q > 0:5, 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 q > 0:5. It is concluded that mixing input power could be halved (viz., q ¼ 0:5) without significant alterations of mixing efficiency. The power input for a single nozzle can be expressed as (Wu ): which corresponds to 1.079 W m À3 ≃ 1 W m À3 for q ¼ 0:5. 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 q ¼ 0:2, but is clearly present in both q ¼ 0:5 and q ¼ 1. Thus, organic material is able to reach the bottom of the tank in both q ¼ 0:5 and q ¼ 1, but not in q ¼ 0:2. This means that: (i) q ¼ 0:2 should be discarded; and (ii) q ¼ 0:5 and q ¼ 1 produce similar mixing effects, and hence, q ¼ 0:5 should be preferred as it incurs less energy consumption.
Optimization strategy second step: finding the best mixing system Throughout the following part of this work, the value of q ¼ 0:5 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 undertake 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 b).
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 nonzero 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 uniformity index texts for both the seedings. Both the seedings show a clear, advantage of the switched strategies over the nonswitched for all the rheology models. Also, an advantage of the '1 min' strategy over the '5 min' is clearly identifiable. Figure 10 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 ), it was claimed in more recent research (Hurtado et al. ) 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.

Flowchart for a general strategy
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 12 referred to this general case as 'altered geometry/injection'.

CONCLUSIONS
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 Uniformity Index 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