Partial integration of ADM1 into CFD: understanding the impact of diffusion on anaerobic digestion mixing

Generally, it is assumed that mixing enhances the performance of anaerobic digestion (AD), creating a homogeneous distribution of soluble substrates, biomass, pH, and temperature. The reported computational fluid dynamics (CFD) models of mechanically agitated anaerobic digesters mostly dealt solely with hydrodynamics. They focused on understanding the effects of sludge rheology, impeller speed, turbulence model selection, geometric configuration, analysis of flow field and evaluation of dead volume for scale-down, laboratory, and pilot-scale anaerobic digesters (Wu & Chen 2008; Meroney & Colorado 2009; Wu 2010; Yu et al. 2011; Bridgeman 2012; Craig et al. 2013; Leonzio 2018; Wiedemann et al. 2018; Conti et al. 2019).

After the development of Anaerobic Digestion Model no. 1 (ADM1) by the International Water Association (IWA) task group, several slight modifications have been formulated and applied in kinetics modelling while little effort has been made to integrate the kinetic model into the CFD model. For example, Wu (2012) modelled an integrated mixing, heat transfer, and fermentation of an egg-shaped anaerobic digester, and showed the distribution of biomass, pH, temperature, mixing, and heat transfer affect methane yield. Gaden & Bibeau (2013) and Rezavand et al. (2019) suggested the implementation of ADM1 into the CFD model, but the stiffness of the kinetic model and solver efficiency proved to be challenging to apply the CFD-kinetic modelling extensively. Other authors like Donoso-Bravo et al. (2018) proposed the integration of ADM1 into a compartmental model intending to understand the spatial variation of soluble substrates, biomass, and pH.

The findings of the AD CFD model hydrodynamics in the literature show that the velocity distribution decreases with increasing total suspended solids (TSS), hence reducing system performance. However, data collected from 2015 to 2017 from a full-scale anaerobic digester in Breda, The Netherlands, indicate that the digester was normally operating at sludge TSS ranging from 3 to 10% and at constant mixing speed. These two contradicting cases suggest that there are knowledge gaps in currently used models that study AD mixing. We hypothesize that the CFD hydrodynamic models alone cannot predict the mixing of AD accurately since the effects of (molecular and turbulent) diffusion and biogas bubbles are ignored. Moreover, velocity is not a good descriptor of mixing homogeneity since minimum velocity required to produce uniform mixing is still unknown. Understanding the effect of diffusion-induced transport on AD mixing and describing the homogeneity of AD in terms of concentration distribution are hence the aim of this work.

This work has two sections. The first section deals with the CFD hydrodynamics modelling of AD. In this section, the homogeneity of AD mixing is described based on velocity distribution. Subsequently, user defined scalar (UDS) transport equations are implemented over the converged and frozen hydrodynamic model in the second section. Part of the ADM1 model is implemented as UDS to simulate concentrations and to allow the effect of diffusion-induced transport to be included. We restricted this to hydrolysis of carbohydrate, protein, and lipids to soluble sugar, amino acids, and long-chain fatty acids (LCFAs), respectively, and growth of biomass, which degrades the hydrolysis products.

The anaerobic digester under study is located in Breda, The Netherlands, and has a volume of 9,000 m³, two stage mixers with three hydrofoil type impellers on each stage, and operates at TSS of 3–10%. The digester was scaled down by the ratio of 1–20 to reduce the computational burden. Figure 1 (left) shows the scaled-down geometry and its detailed dimensions.

A multiple reference frame is selected to implement the flow feature of the impeller rotation. A rotating reference frame (RRF) for each stage of the impellers and a stationary reference frame are created. The impellers are enclosed in an RRF and connected to the stationary reference frame using an interface boundary condition.

The stationary and rotating reference frames are meshed by a sweeping method for sweepable bodies and hex-dominant for non-sweepable bodies. The mesh independence test runs for three different mesh sizes consisting of 1.2, 1.6, and 3.2 million elements. Figure 1 (right) shows the final mesh with 3.2 million cells, for which the simulation proved to be independent of mesh size.

In this work, 4% of TSS sludge is considered to study the impact of diffusion-induced transport. The occurrence of density gradients due to digested sludge sedimentation is expected to be very low. For this reason, as well as to avoid the complexity of a multiphase model introducing many additional empirical parameters to describe the momentum transfer between the phases, simulating the system as a single-phase fluid is reasonable. The free surface at the top of the digester is ignored and treated as a wall boundary condition. An inlet velocity is defined at the inlet of the digester (Figure 1). Since the details of the flow velocity and pressure are not known before solving the flow problem, the outflow boundary condition is defined at the outlet. The hydrodynamics of AD mixing is solved, treating the sludge rheology as non-Newtonian fluid, employing a Herschel–Bulkley rheology model and k- turbulence model. The rotational speed of both mixers is set to 50 rpm.

The scalar variable (⁠ in Equations (1) and (2) results from either hydrolysis (Equation (7)) or growth of biomass on soluble substrates (Equation (8)). Three hydrolysis products, three biomass growth, and two diffusion equations, as well as three initial conditions are written in C and compiled into Ansys Fluent 19 R1.

Since the sludge compressibility is very low, the pressure-based solver SIMPLE (Semi-Implicit Method for Pressure Linked Equations) is selected. Momentum, turbulent kinetic energy, turbulence dissipation, and all UDS scalar variables are discretized using a first-order upwind method while the pressure is discretized using the standard method. A computer cluster of 4 × AMD Opteron 6380 2.5 GHz (16 cores) is used for the simulation.

Figure 2(a) shows the contour plot of velocity along the vertical plane. The figure indicates that only sludge near the impeller is moving at significant velocity, whereas inlet velocity does not produce significant velocity compared to the velocity produced by the impeller, and hence its impact on AD mixing is insignificant. The velocity produced by the top and bottom impellers shows a significant difference even though both impellers rotate at the same mixing speed. The reason behind this is explained next, together with the velocity vector shown in Figure 2(b).

Figure 2(b) shows that the velocity vectors produced by the top and bottom impellers point in opposite directions. Hence, the sum of the top and bottom impellers' velocity produces a relative velocity which cancels each other, defeating the purpose of having two impellers, and resulting in low-velocity distribution in the vicinity of the top impeller. Furthermore, this makes a design based on two axial flow impellers in the same axis questionable, and a combination of an axial and radial impeller would make more sense.

In addition to the velocity contour and velocity vector plot, non-dimensional velocity u to the maximum velocity at the tip of the impeller, ⁠, i.e. along r/R at different height of the digester, is plotted in Figure 2(c). The figure indicates approaches 1 near the tip of the impellers and approaches 0 near the center and wall of the digester. The variation of indicates that velocity dissipates quickly due to high sludge viscosity as it moves away from the impeller tip.

The cumulative volume distribution of velocity elegantly summarizes the 3-D velocity distribution (Figure 2(d)). The figure quantifies that velocity distribution across the digester volume is not uniform. For example, about 75% of digester volume has a velocity of <5% of maximum velocity, which indicates that the velocity variation between the tip of the impeller and the rest of the digester volume is very high.

Mesh independence of the simulation is analyzed based on the mass balance between the inlet and outlet of the digester and the evaluation of cumulative velocity distribution variation with mesh size (Figure 2(d)). A closed mass balance is attained in all cases, and variation of velocity distribution in an anaerobic digester volume is not significant as such. For example, the cumulative velocity distribution variation between 1.6 million and 3.2 million elements is less than 7%.

Figures 3 and 4 show contour plots of soluble sugar concentration and biomass distribution variation under advection and different advection–diffusion transport. This is further clarified by the Pe variation tabulated in Table 1. The table indicates that Pe varies along r/R and height (bottom to top) of the digester. Regions around the top, bottom, near the corners, and walls of the digester have a Pe < 1, which means diffusion transport is dominant in these regions. The remaining regions have a Pe > 1 meaning that advection transport is dominant.

Figures 3(a) and 4(a) indicate the contour plots of soluble sugar and biomass concentration under advection transport (Equation (2)). Here, the concentration of soluble sugar and biomass varies widely between the regions with Pe > 1 and Pe < 1. The maximum concentrations of soluble sugar and biomass are 39 g/m³ and 3.6 g/m³, respectively, in a region with Pe < 1, and the corresponding minimum concentrations are 30 g/m³ and 2.4 g/m³, respectively, in a region with Pe > 1. Due to the high concentration of soluble sugar and biomass near the corners and walls, the soluble sugar and biomass concentration distributions in a region with Pe > 1 are low, hence indicating poor mixing uniformity.

Figures 3(b) and 4(b) show the distribution of soluble sugar and biomass under advection–molecular diffusion transport, respectively. In this case, the homogeneity of soluble sugar and biomass concentration in a region with Pe > 1 increases because of molecular diffusion from a region with Pe < 1 to a region with Pe > 1. The lowest soluble sugar and biomass concentration, distribution along the plane in Figures 3(b) and 4(b) changed to 33 g/m³ and 2.7 g/m³, respectively, under advection–molecular diffusion transport, corresponding to 30 g/m³ and 2.7 g/m³ under sole advection transport. The highest concentration near the corner with Pe < 1 is still 39 g/m³ and 3.6 g/m³ under advection–molecular diffusion transport.

The impact of mixing due to advection–turbulence diffusion is shown in Figure 3(c)–3(e) for different Sc_t. Under advection–turbulent diffusion, the uniformity of soluble sugar concentration is much higher compared to the advection and advection–molecular diffusion transport. The concentration distribution for Sc_t = 0.9 and 0.7 is almost similar. The lowest concentration of soluble sugar is about 35 g/m³ for Sc_t = 0.9 and 0.7, while it is about 37 g/m³ for Sc_t = 0.5. Uniformity of soluble sugar increases with decreasing Sc_t, and this indicates that turbulent mixing efficiency increases when turbulent eddy diffusivity is dominating turbulent viscosity (Equation (5)).

In the advection–(molecular and turbulent) diffusion transport model (Figures 3(f) and 4(c)), the homogeneity of both soluble sugar and biomass concentrations was found to be higher than all cases discussed. Unlike other models, in advection–(molecular and turbulent) diffusion transport, the ranges of minimum and maximum concentration distribution in a region with Pe > 1 and Pe < 1 are much more comparable. For example, the minimum and maximum concentrations of soluble sugar are about 37 g/m³ and 39 g/m³, respectively. Since both molecular and turbulent diffusion always exist together, it is recommended to model combined advection–diffusion (molecular and turbulent) to describe AD mixing homogeneity accurately. The contour plot and cumulative volume distribution discussed did not include the results of soluble amino acids, LCFAs, and their corresponding biomass because their concentration distribution patterns are similar to soluble sugar and biomass.

In a similar approach to non-dimensional velocity distribution shown in Figure 2(b), the non-dimensional concentration distribution of soluble sugar is plotted along r/R at different heights of the digester for advection and advection–diffusion transport models (Figure 5). Soluble sugar concentration is non-dimensionalized by dividing the concentration of soluble sugar along r/R by the maximum concentration of the soluble sugar (⁠⁠) at the wall. The ratio of and (⁠ along r/R at different height of the digester ranges between 0.6 and 1 (⁠1) for the advection transport model (Figure 5 left).

Figure 5 (right) shows is between 0.96 and 1 (⁠1) at all indicated heights, indicating the digester is nearly ideally mixed under advection–diffusion transport. Unlike non-dimensional velocity distribution, (Figure 2(c)), non-dimensional concentration distribution, (Figure 5), along r/R is uniform. So, making conclusions about AD mixing homogeneity based on velocity and concentration distribution leads to very different outcomes.

In addition to concentration distribution contour plots and a non-dimensional plot along r/R and height (Figures 2–5) in 2-D, additional information can be extracted by plotting a concentration distribution derived from the 3-D volume.

A comparison is provided for soluble substrates and biomass sole advection transport as well as all advection–diffusion transport cases modelled. Distributions of soluble sugar and biomass shift to the right for all advection–diffusion transport cases compared to the advection-only transport model. The shaded area under the curve indicates the increase in the uniformity of AD mixing by advection–diffusion for that case. For example, the shaded area in Figures 6(a) and 7(a) shows a change of soluble sugar and biomass concentration distribution, respectively, for advection–molecular diffusion, and the shaded area is small. A small shaded area means the difference of soluble substrates and biomass distribution under advection–molecular and sole advection transport is rather small. Mainly, it is minimal for biomass distribution due to the small molecular diffusion coefficient of biomass (Figure 7(a)).

Figure 6(b)–6(d) show the effect of advection–turbulent diffusion transport for different Sc_t. The figures indicate that the shaded area increases with decreasing the Sc_t, which means mixing performance increases when eddy mass diffusivity dominates turbulent diffusion (Equation (5)). The shaded area under advection–turbulent diffusion transport is larger than under advection–molecular diffusion transport. The shaded area is almost the same for a model with Sc_t = 0.7 and Sc_t = 0.9, and the largest for a model with Sc_t = 0.5.

In advection–(molecular and turbulent) diffusion transport, the homogeneity of mixing increased much more compared to the rest of the advection–diffusion transport cases considered (Figures 6(e) and 7(b)) except advection–turbulent diffusion with Sc_t = 0.5. The shaded area comparison of advection–molecular diffusion, advection–turbulent diffusion, and advection–(molecular and turbulent) diffusion transport shows that advection–diffusion yields a higher homogeneity of AD.

Generally, the steeper the slope of the cumulative volume distribution, the better the homogeneity of AD. In other words, change in concentration distribution variation is smaller, with change in cumulative volume distribution at the steepest slope.

In the conventional CFD hydrodynamics model, the uniformity of mixing comparison is based on relative velocity distribution (Figure 2), i.e. a high-velocity region near the impeller tip is taken as a well-mixed region, which is used as a reference point to compare the remaining velocity distribution in an anaerobic digester. The analysis of velocity distribution leads to the conclusion that the lowest velocity region relative to the velocity near the impeller tip is assumed as not mixed well irrespective of mixing due to turbulence and molecular diffusion. So, making conclusions about the homogeneity of AD in terms of relative velocity distribution analysis is tricky and can lead to inaccurate conclusions. Because the minimum velocity required to produce homogeneity is not known and the conventional CFD model ignores the contribution of (molecular and turbulence) diffusion.

In the advanced advection–diffusion transport model, the limitations of the conventional CFD model mixing description improved in two ways. First, the homogeneity of mixing is described based on concentration distribution. Second, the effects of mixing due to diffusion transport are included. Describing the homogeneity of AD based on concentration distribution, including the effect of diffusion transport on AD mixing, gives more comprehensive information.

The results of different advection–diffusion transport model cases shown in the contour plot (Figures 3 and 4) and cumulative volume distribution (Figures 6 and 7) shows that changing the diffusion transport variables affects the homogeneity of soluble substrates and biomass. The models indicate that soluble substrates and biomass are less homogeneous under the sole advection transport, and the homogeneity of soluble substrates and biomass increases under advection–diffusion transport.

The cumulative volume distribution under advection–diffusion transport (Figures 6(f) and 7(b)) indicate that the variation of concentration distribution is minimal, and it is close to a homogeneously mixed AD. This shows that the effect of diffusion transport is significant in the homogenization of AD and should be considered in AD mixing optimization. Since molecular diffusion is a material property and mixing optimization cannot improve it, considering turbulence generation techniques in anaerobic digester and mixer design increases the homogeneity of substrates and biomass distribution.

In general, CFD hydrodynamics and advection–diffusion transport give significantly different results. Velocity is a carrier/transporter of scalar variables like substrates and biomass, and it is not a good mixing descriptor by itself. So, explaining the homogeneity of AD based on concentration improves the interpretation of AD mixing, including the impact of diffusion transport.

Conventional CFD hydrodynamics models solely based on the velocity description fall short of describing the homogeneity of AD mixing. Advanced CFD modelling, including advection–diffusion transport, improves AD mixing description and understanding. Mixing seems to be more profound based on concentration profiles than what would be expected from velocity distributions. Analyzing the homogeneity of AD mixing in terms of concentration significantly facilitates the interpretation of the computed results. As the contribution of (molecular and turbulent) diffusion transport is significant, it should be included in future mixing optimization studies.

This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 676070.

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