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

Suf ﬁ cient mixing is crucial for the proper performance of anaerobic digesters (ADs), creating a homogeneous distribution of soluble substrates, biomass, pH, and temperature. The opaqueness of the sludge and mode of operation make it challenging to study AD mixing experimentally. Therefore, hydrodynamics modelling employing Computational Fluid Dynamics (CFD) is often used to investigate this mixing. However, CFD models mostly do not include biochemical reactions and, hence, ignore the effect of diffusion-induced transport on AD heterogeneity. The novelty of this work is the partial integration of Anaerobic Digestion Model no. 1 (ADM1) into the CFD model. The aim is to better understand the effect of advection-diffusion transport on the homogenization of soluble substrates and biomass. Furthermore, AD homogeneity analysis in terms of concentration distribution is proposed rather than the traditional velocity distributions. The computed results indicate that including diffusion-induced transport affects the homogeneity of AD.


INTRODUCTION
Generally, it is assumed that mixing enhances the performance of AD, creating a homogeneous distribution of soluble substrates, biomass, pH, and temperature. The reported CFD models of mechanically agitated ADs 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 scaledown, laboratory, and pilot-scale AD (Conti et  After the development of ADM1 by the International Water Association (IWA) task group, several slight modifications have been formulated and applied in kinetics modelling while little effort have been made to integrate the kinetic model into CFD model. For example, Wu () modelled an integrated mixing, heat transfer, and fermentation of an egg-shaped AD, and showed the distribution of biomass, pH, temperature, mixing, and heat transfer affect methane yield. Gaden & Bibeau () and Rezavand et al. () 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. () 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 AD in Breda, the Netherlands, indicates 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.

METHODS
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 allowing to include the effect of diffusion-induced transport. We restricted this to hydrolysis of carbohydrate, protein, and lipids to soluble sugar, amino acids, and Long Chain-Fatty Acid (LCFA), respectively, and growth of biomass, which degrades the hydrolysis products are part of ADM1 integrated into CFD model.

GEOMETRY, MESH AND HYDRODYNAMICS MODELLING
The AD under study is located in Breda, The Netherlands, and has a volume of 9,000 m 3 , 2 stage mixers with three hydrofoil type impellers on each stage, and operates at TSS 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 (MRF) 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 RRF are meshed by a sweeping method for sweepable bodies and Hex-dominant for nonsweepable 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 total suspended solids (TSS) sludge is considered to study the impact of diffusioninduced 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 singlephase 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.

GOVERNING EQUATIONS OF UDS TRANSPORT MODELLING
The UDS transport equations are solved after the hydrodynamics simulation has converged. Continuity, momentum, and turbulence modelling equations are disabled (frozen) before implementing and solving the UDS transport equations. Activating the source term in ANSYS FLUENT enables advection-diffusion Equation (1) with diffusion transport or Equation (2) without diffusion transport for steady state simulation.
where ϕ k is the scalar variable, k ¼ 1, 2,…, N. S ϕ k is the source term, ρ is the sludge density and Γ eff is the effective (molecular and turbulent) diffusion coefficient (Equation (3)).
where D i,lam is the molecular diffusion coefficient. Due to a lack of substrates and biomass diffusion coefficients for sludge, molecular diffusion coefficients (D i,lam ) of glucose in an alginate gel solution 6.4*10 À10 m 2 /s are adopted for all substrates (Øyaas et al. ). Furthermore, a bacteria diffusion coefficient 2.8*10 À11 m 2 /s is adopted for all biomass (Douarche et al. ) and μ trepresents the turbulent viscosity (Equation (4)).
where k is the turbulent kinetic energy, ε is the eddy dissipation and the details of C μ is found in Ansys (). Sc t is the turbulent Schmidt number, which is an empirical constant and relatively insensitive to the fluid properties. Sc t is the ratio of turbulent momentum diffusivity (eddy viscosity), v t to mass diffusivity, D t (Equation (5).
Usually, the value of Sc t is between 0.5 and 0.9 (Tominaga & Stathopoulos ). The effect of turbulent diffusion on AD mixing is analyzed for Sc t values, i.e. 0.5, 0.7, and 0.9.
Introducing the non-dimensionless Peclet (Pe) number (Equation (6)) helps to describe advection-diffusion transport clearly. Pe is the ratio of advection transport to diffusion transport. In regions where Pe < 1, diffusion transport is dominant while in regions with Pe > 1, advection transport is dominant (Ghatak ).
where u is the flow velocity, d is the characteristic length and D is the molecular diffusion coefficient.

ADM1 EQUATIONS IMPLEMENTED INTO CFD MODEL
The scalar variable (ϕ k ) 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. The first-order rate Equation (7) represents the hydrolysis of carbohydrates, lipids, and proteins to their respective soluble substrate.
where S j is the product of hydrolysis, i.e., soluble sugars, amino acids and LCFAs from X i , i.e. carbohydrates, proteins, and lipids, respectively. The Monod Equation (8) is used to represent the growth of biomass on each of these soluble substrates.
where Y j is the yield of biomass on a substrate, I pH is the pH inhibition (0.992) and I IN,lim is the nitrogen limitation inhibition (0.998) and k m,j is the maximum specific growth rate. The values of biomass yield from the substrates (Y j ), the hydrolysis constants (k hyd,i ), half-saturation value (K S,j ), and Monod maximum specific growth rate (k m,j ) are taken from Rosen & Jeppsson () and are tabulated in Table A1 in the appendix. For further details on ADM1 and its implementation the reader is referred to Batstone et al. () and Rosen & Jeppsson (), respectively. The initial concentration of carbohydrates, lipids, and proteins are taken as 19.5 g/L, 45.1 g/l, and 8.1 g/L, respectively. 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 1 st order upwind method while the pressure is discretized using the standard method. A Computer Cluster of 4 x AMD Opteron 6380 2.5 GHz (16 Cores) is used for the simulation.

CFD hydrodynamics model
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   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, U tip , i.e., u U tip along r/R at different height of the digester, is plotted in Figure 2(c). The figure indicates u U tip approaches one near the tip of the impellers and approaches zero near the center and wall of the digester. The variation of u U tip 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 AD 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%.
Advanced CFD-kinetics modelling: a comparison of advection and advection-diffusion transport   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 concentration of soluble sugar and biomass are !39 g/m 3 and !3.6 g/m 3 , respectively, in a region with Pe < 1, and the corresponding minimum concentrations are 30 g/m 3 and 2.4 g/m 3 , 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 distribution 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 3 and 2.7 g/m 3 , respectively, under advection-molecular diffusion transport corresponding 30 g/m 3 and 2.7 g/m 3 under sole advection-transport. The highest concentration near the corner with Pe < 1 is still !39 g/m 3 and !3.6 g/ m 3 under advection-molecular diffusion transport.

Concentration distribution contour plot
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 3 for Sc t ¼ 0.9 and 0.7 while it is about 37 g/m 3 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 concentration were found higher than all cases discussed. Unlike other models, in advection-(molecular and turbulent) diffusion transport, the range 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 concentration of soluble sugar are about 37 g/m 3 and 39 g/m 3 , respectively. Since both molecular and turbulent diffusion always exists together, it is recommended to model combined advection-diffusion (molecular and turbulent) to describe AD mixing homogeneity accurately.

Non-dimensional concentration variation along the radius and height
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 concentration (C) along r/R by the maximum concentration of the soluble sugar concentration (C wall ) at the wall. The ratio of C and C wall ( C C wall ) along r/R at different height of the digester ranges between 0.6 and 1 (0:6 C C wall 1) for the advection transport model (Figure 5 left).

Mixing uniformity analysis using cumulative volume distribution
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 and this for 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. 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) shows 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 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 advectionturbulent 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 steepest the slope of the cumulative volume distribution, the better the homogeneity of AD. In other words, change in concentration distribution variation

DISCUSSION
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, and it uses as a reference point to compare the remaining velocity distribution in AD. 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 homogenous under the sole advection transport, and the homogeneity of soluble substrates and biomass increases under advection-diffusion transport.
The cumulative volume distribution under advectiondiffusion 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 AD 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.

CONCLUSION
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.