Linking CFD and Kinetic Models in Anaerobic Digestion Using a Compartmental Model Approach

: Understanding mixing behavior and its impact on conversion processes is essential for the operational stability and conversion e ﬃ ciency of anaerobic digestion (AD). Mathematical modelling is a powerful tool to achieve this. Direct linkage of Computational Fluid Dynamics (CFD) and the kinetic model is, however, computationally expensive, given the sti ﬀ ness of the kinetic model. Therefore, this paper proposes a compartmental model (CM) approach, which is derived from a converged CFD solution to understand the performance of AD under non-ideal mixing conditions and with spatial variation of substrates, biomass, pH, and speciﬁc biogas and methane production. To quantify the e ﬀ ect of non-uniformity on the reactor performance, the CM implements the Anaerobic Digestion Model 1 (ADM1) in each compartment. It is demonstrated that the performance and spatial variation of the biochemical process in a CM are signiﬁcantly di ﬀ erent from a continuously stirred tank reactor (CSTR) assumption. Hence, the assumption of complete mixed conditions needs attention concerning the AD performance prediction and biochemical


Introduction
Anaerobic digestion (AD) is a widely used technology in wastewater sludge management, because of its ability to convert organic compounds in the sludge to renewable energy (biogas) to reduce sludge volume and to remove odors related to the sludge. In practice, mixing is conducted as it is believed to improve the process performance by creating a homogeneous distribution of soluble substrates, biomass, pH, and temperature.
The AD modelling in the literature can be categorized into two categories. The first model type is kinetic modelling of AD which assumes the digester is a continuously stirred tank reactor (CSTR) and tries to understand the AD kinetics, ignoring the hydrodynamics and sludge properties, by using the Anaerobic Digestion Model 1 (ADM1), developed by International Water Association (IWA) specialist group [1].
Several kinetic models of AD have been reported based on the ADM1 model. For example, Polizzi et al. [2] simulated primary tannery sludge; Kil et al. [3] simulated the relationship between the steady-state process and concentration of the influent components to estimate the output waste characteristics of wastewater sludge. Poggio et al. [4] used the ADM1 to simulate the digestion of solid organic waste, and Donoso-Bravo et al. [5] implemented the ADM1 into a compartmental model to understand the spatial variation of soluble substrates, biomass, and pH in AD based on a simple division of the AD volume into sub-volumes. The MRF has two flow features: a rotating reference frame (RRF), which encloses the volume surrounding the mixer impellers, and the stationary reference frame. An RRF was connected to the stationary reference frame using an interface boundary condition. The impeller speeds used in the ADs range from low speed to very high speed (up to 800 rpm) depending on the sludge TSS [9,15,24]. We selected an intermediate value of 300 rpm.
Both the stationary and rotating frames of reference were meshed using a tetrahedral mesh with 1.66 million elements (cells) (0.7 mm 3 /cell) and 0.312 million nodes (Figure 1, left) in which the results are independent of mesh sizes at residual of 10 −3 . A mesh independent test was analyzed for 95,000 cells, 1.3 million cells, and 1.66 million cells. The percentage error of the velocity distribution between 1.35 and 1.66 million cells is about 5%. Inlet velocity, pressure outlet, and no-slip wall were defined for inlet, outlet, and wall boundary conditions, respectively.

Discretization Method
The SIMPLE (Semi-Implicit Method for Pressure Linked Equations) scheme solver was used for the solution of a steady-state with a k-ε realizable turbulence model. The momentum and turbulent kinetic energy were discretized using a 2nd order upwind scheme while the turbulence dissipation was discretized using a 1st order upwind scheme since it oscillates and does not converge when using a 2nd order upwind scheme.
Continuity and momentum equations of an RRF for absolute velocity formulation are given in Equations (1) and (2), respectively.
where → is the relative velocity (the velocity observed from the rotating frame (m/s)), ω → is the angular velocity (1/s), and → is the absolute velocity viewed from the stationary reference frame (m/s).
The realizable turbulence kinetic energy (k) and dissipation rate (ε) were selected based on the work of Wu [25] and given by Equations (3) and (4), respectively.

Discretization Method
The SIMPLE (Semi-Implicit Method for Pressure Linked Equations) scheme solver was used for the solution of a steady-state with a k-ε realizable turbulence model. The momentum and turbulent kinetic energy were discretized using a 2nd order upwind scheme while the turbulence dissipation was discretized using a 1st order upwind scheme since it oscillates and does not converge when using a 2nd order upwind scheme.
Continuity and momentum equations of an RRF for absolute velocity formulation are given in Equations (1) and (2), respectively.
where → v r is the relative velocity (the velocity observed from the rotating frame (m/s)), → ω is the angular velocity (1/s), and → v is the absolute velocity viewed from the stationary reference frame (m/s). The realizable turbulence kinetic energy (k) and dissipation rate (ε) were selected based on the work of Wu [25] and given by Equations (3) and (4), respectively.
where C 1 = max 0.43, η η+5 , η = S k ε , S = 2S ij S ij , C 1ε , C 2 , σ k , σ ε are constants with a value of 1.44, 1.9, 1.0 and 1.2, respectively, G k is the generation of turbulent kinetic energy due to mean velocity gradients, G b is the generation of turbulent kinetic energy due to buoyancy, S is the modulus of the strain rate tensor, and µ t is eddy viscosity. The detailed explanation of Y M , G k , S and C µ are found in the ANSYS Theory Guide [26].
The Herschel-Bulkley rheology model was selected because it includes non-zero shear stress (yield stress) even when the strain rate is zero [27]. Hence, Equation (5)  γ c . A non-Newtonian rheology model was activated by the user using the Text User Interface (TUI) in the ANSYS FLUENT console.
where η is the apparent viscosity (Pa. s), .
γ is the shear rate (1/s), K is the consistency index, n is the flow index, and . γ c is the critical shear rate. The critical shear rate was added to avoid discontinuity (infinite apparent viscosity) at a very low shear rate [26,27]. K is the consistency index (Pa.sn), n is the power-law index and τ y is the yield stress (Pa). The values of K, n and τ y were correlated based on the work of Eshtiaghi et al. [28] and Slatter [29] and are shown in Table 1. The effect of sludge rheology was modelled by treating the sludge as a Newtonian fluid (water) for sludge with TSS less than 4% and non-Newtonian fluid with the TSS at 4%, 8%, and 12%.

Compartmentalization of AD from the CFD Model
A CM is a general method of dividing a reactor volume (here an AD) into sub-volumes based on the CFD model velocity distribution to mimic the non-ideal mixing kinetics model. Different methods of compartmentalization have been discussed in the literature, depending on the type of problem. For example, velocity, species concentration, turbulence, and temperature can be applied to compartmentalize the system under study [22] using either manual [20,30] or automatic zoning [31] of the system under study. Due to its flexibility, simplicity, and the fact that the user can decide the number of compartments, the digester was compartmentalized manually based on the velocity distribution in this work. Figure 2 shows a method of deriving the CM from a CFD hydrodynamic model and links it to AD kinetics. The compartments are linked together by the flux exchange arrows.
The exchange flux between each compartment depends on the local velocity magnitude and the compartment surface area (Equation (7)).
where u avg is the average velocity (m/s) in/out of the compartment and A sur f is the surface area of the compartment (m 2 ). In a CM, each compartment is assumed to be ideally mixed, and as having virtually independent digesters; hence, the ADM1 is implemented in each compartment. The exchange flux between each compartment depends on the local velocity magnitude and the compartment surface area (Equation (7)).
where is the average velocity (m/s) in/out of the compartment and A is the surface area of the compartment (m 2 ).
In a CM, each compartment is assumed to be ideally mixed, and as having virtually independent digesters; hence, the ADM1 is implemented in each compartment.

Implementation of ADM1 into the CM
The steady-state mass balance of the ADM1 soluble substrates and biomass are given in Equations (8) and (9) [1].
where q and q are the inflow and outflow rate, respectively, S represents the soluble substrates, X represents the biomass fractions, V is the digester liquid volume, is the uptake rate of process , v , represents the rate coefficients for component i taking part in process , represents the index soluble substrates and biomass index, and represents the number inflow rates into the compartment. Since there might be more than one inflow rate into a compartment, the total inflow rate is indicated as the sum of inflow rates.
The biogas flow rate ( ) from each compartment at atmospheric pressure was calculated using Equation (10).
where is related to the biogas outlet friction factor with 5×10 4 m 3 d −1 bar −1 and is the sum of the partial pressure of methane, CO2, and water vapor, and Patm is the atmospheric pressure.
The CM shown in Figure 3 was implemented in WEST software (powered by DHI, Denmark) after extensive customization of the software to fit the modelling goal. The ADM1 model parameter

Implementation of ADM1 into the CM
The steady-state mass balance of the ADM1 soluble substrates and biomass are given in Equations (8) and (9) where q in and q out are the inflow and outflow rate, respectively, S liq represents the soluble substrates, X liq represents the biomass fractions, V liq is the digester liquid volume, ρ j is the uptake rate of process j, v i,j represents the rate coefficients for component i taking part in process j, i represents the index soluble substrates and biomass index, and k represents the number inflow rates into the compartment. Since there might be more than one inflow rate into a compartment, the total inflow rate is indicated as the sum of inflow rates. The biogas flow rate (q gas ) from each compartment at atmospheric pressure was calculated using Equation (10). q gas = k p P gas − P atm P gas P atm (10) where k p is related to the biogas outlet friction factor with 5 × 104 m 3 d − 1 bar − 1 and P gas is the sum of the partial pressure of methane, CO 2, and water vapor, and P atm is the atmospheric pressure. The CM shown in Figure 3 was implemented in WEST software (powered by DHI, Denmark) after extensive customization of the software to fit the modelling goal. The ADM1 model parameter values were taken as proposed in the Benchmark Simulation Model 2 (BSM2) framework [32]. The ADM1 input data were fractionated from chemical oxygen demand (COD), Total kjeldahl Nitrogen (TKN, total suspended solids (TSS), and inflow rate, according to Nopens et al. [33]. The input values of COD, TKN, TSS, and inflow rate were 66 g/L, 0.8 g/L, 47 g/L, and 75 L/d, respectively. uniformity decreases sharply for non-Newtonian fluid as the sludge TSS increases. At an elevated TSS, only the sludge near the impeller is moving due to higher sludge viscosity. This results in the poor velocity distribution ( Figure 3) in a similar fashion reported by [9].
Since the contour plot is a qualitative observation and only gives observation in a plane view in 2-D, quantitative velocity distribution in 3-D is necessary for detailed information on the uniformity of velocity distribution. Quantitative velocity distribution in 3-D is plotted using the cumulative volume frequency of velocity shown in Figure 3 (b) for Newtonian and non-Newtonian fluid.

CFD Model Simulation
The steady-state simulation results of the AD mixing treating the sludge as Newtonian and non-Newtonian fluid are shown in Figure 3a in a vertical plane view. The observation from Figure 3a indicates that the velocity distribution has a good uniformity for Newtonian fluid (water), and the uniformity decreases sharply for non-Newtonian fluid as the sludge TSS increases. At an elevated TSS value, only the sludge near the impeller is moving due to a higher sludge viscosity. This results in a poor velocity distribution ( Figure 3) in a similar fashion reported by [9].
Since the contour plot is a qualitative observation and only allows for observation in a plane view in 2-D, quantitative velocity distribution in 3-D is necessary for detailed information on the uniformity of velocity distribution. Quantitative velocity distribution in 3-D is plotted using the cumulative volume frequency of velocity shown in Figure 3b for Newtonian and non-Newtonian fluid.
The maximum velocity near the impeller speed is 2.3 m/s; however, the volume of the sludge, which moves above 1 m/s is minimal and cannot be observed in the cumulative volume frequency of velocity distribution analysis and contour plot. Based on the assumption stated in [9], the digester zone moving at a velocity of less than 5% of the maximum velocity is considered as a stagnant zone. Hence, 10%, 58%, 93%, and 97% of the digester volume has a velocity of less than 5% of the maximum velocity for water and sludge with the TSS at 4%, 8%, and 12%, respectively.
In addition to the vertical plane velocity contour plot and cumulative volume frequency of velocity, the velocity distribution in a horizontal plane view is plotted at different heights of the digester for the sludge with the TSS at 4% (Figure 3d) at an indicated height in Figure 3c to obtain a precise observation  Figure 3c indicates a velocity contour plot of the sludge with the TSS at 4% and sliced at different heights. Since the impellers are axial flow types, the sludge flows through the impeller, and downward velocity is created as indicated in the velocity vector plot (Figure 4). Then the velocity bounces at the bottom of the digester and creates two 'U' shaped patterns (Figure 3c), and the uniform velocity distribution is observed at h 1 (Figure 3d). The magnitude and uniformity of velocity decrease from h 1 to h 5 . There is a high velocity around the impellers, medium velocity near the wall, and low velocity between the wall and center of the digester. The reason for such a flow formation is due to the nature of the impeller. From h 5 to the top of the digester, the velocity distribution gradually decreases and becomes stagnant. In addition to the vertical plane velocity contour plot and cumulative volume frequency of velocity, the velocity distribution in a horizontal plane view is plotted at different heights of the digester for the sludge with the TSS at 4% (Figure 3d) at an indicated height in Figure 3c to obtain a precise observation of the velocity distribution together with the contour plot in the vertical plane view (Figure 3a) in zoning the digester. Figure 3c indicates a velocity contour plot of the sludge with the TSS at 4% and sliced at different heights. Since the impellers are axial flow types, the sludge flows through the impeller, and downward velocity is created as indicated in the velocity vector plot (Figure 4). Then the velocity bounces at the bottom of the digester and creates two 'U' shaped patterns (Figure 3c), and the uniform velocity distribution is observed at h1 (Figure 3d). The magnitude and uniformity of velocity decrease from h1 to h5. There is a high velocity around the impellers, medium velocity near the wall, and low velocity between the wall and center of the digester. The reason for such a flow formation is due to the nature of the impeller. From h5 to the top of the digester, the velocity distribution gradually decreases and becomes stagnant.
The velocity vector shown in Figure 4 indicates two types of recirculation. Solid circular arrows indicate recirculation between the high and medium velocity zones, and flow exchange between the zone and its neighborhood exists. In contrast, the dashed arrows indicate the recirculation, which creates a closed-loop. The sludge inside this loop recirculates continuously over an extended time because the exchange flow with its neighborhood is limited. A low-velocity recirculation is also observed in the top right side of the stagnant region.

Compartmentalization of AD
The digester is, therefore, categorized into four zones and eight compartments ( Figure 5) similar to the work of Le Moullec [20]. The volume of each compartment is calculated from the cumulative volume distribution of velocity shown in Figure 3b.
1. Zone one: high-velocity zone (C3). This zone is near the impeller with a velocity magnitude ≥ 0.6*maximum velocity (v → ). The local residence time in this zone is very short due to the high flux exchange. v → is the sludge velocity near the mixer impeller tip. 2. Zone two: medium velocity zone. The velocity in this zone ranges between 0.25v This zone is the outermost of the digester volume with a large diameter with low radial velocity; therefore, the recirculation time is longer. Instead of assuming the zone is a single compartment, it is decided to implement two compartments in series (C1 and C5) to make the residence time more realistic. 3. Zone three: low-velocity zone. A zone with a velocity magnitude of 0.05v This zone is divided into two compartments in series (C2 and C4), such as zone two. The velocity vector shown in Figure 4 indicates two types of recirculation. Solid circular arrows indicate recirculation between the high and medium velocity zones, and flow exchange between the zone and its neighborhood exists. In contrast, the dashed arrows indicate the recirculation, which creates a closed-loop. The sludge inside this loop recirculates continuously over an extended time because the exchange flow with its neighborhood is limited. A low-velocity recirculation is also observed in the top right side of the stagnant region.

Compartmentalization of AD
The digester is, therefore, categorized into four zones and eight compartments ( Figure 5) similar to the work of Le Moullec [20]. The volume of each compartment is calculated from the cumulative volume distribution of velocity shown in Figure 3b.

1.
Zone one: high-velocity zone (C3). This zone is near the impeller with a velocity magnitude ≥ 0.6*maximum velocity ( → v max ). The local residence time in this zone is very short due to the high flux exchange.
→ v max is the sludge velocity near the mixer impeller tip.

2.
Zone two: medium velocity zone. The velocity in this zone ranges between 0.25 This zone is the outermost of the digester volume with a large diameter with low radial velocity; therefore, the recirculation time is longer. Instead of assuming the zone is a single compartment, it is decided to implement two compartments in series (C1 and C5) to make the residence time more realistic.

3.
Zone three: low-velocity zone. A zone with a velocity magnitude of 0.05 → v max . This zone is divided into two compartments in series (C2 and C4), such as zone two. 4.
The CM accuracy increases with a higher number of compartments, and it can be derived by reducing the tolerance of velocity variation. Hence, the user should decide the required number of CM in conjunction with the modelling goal. It should be bear in mind that every additional compartment adds 32 differential-algebraic equations (DAEs) to the model.

Non-Ideal Mixing Kinetics Modelling
The spatial variation of biomass, soluble substrates, pH, specific biogas, and methane production in the CM is discussed below and compared with the CSTR AD of the same volume. Figure 6a shows that the biomass concentration in the CM is 35% less than the corresponding CSTR AD. The reduction in hydraulic retention time (HRT) due to the stagnant zone in the CM affects the growth of biomass; hence, it affects the uptake of soluble substrates and reduces the digester performance. The reduction in the biomass concentration in a CM and the stagnant volume (58% of digester volume) are not linearly proportional because a limited amount of biomass grows in the stagnant zones.
The growth of acetogenic biomass (converting volatile fatty acids to acetate) and methanogenic biomass (converting acetate into biogas) are the two most affected in the CM (Figure 6b). The spatial variation of biomass, soluble substrates, pH, specific biogas, and methane production in the CM is discussed below and compared with the CSTR AD of the same volume. Due to the lower biomass growth in the CM (Figure 6a), the undigested soluble substrates outflow from the digester is higher compared to the CSTR digester. The undigested outflow of the soluble acetate and LCFA is 0.4 g/L and 1.4 g/L, respectively, in the CM compared to only 0.04 g/L and 0.1 g/L in the CSTR AD. Moreover, the specific biogas and methane production are reduced by 14% and 16%, respectively, in the CM compared to the CSTR of the same volume (Figure 6c). The CM accuracy increases with a higher number of compartments, and it can be derived by reducing the tolerance of velocity variation. Hence, the user should decide the required number of CM in conjunction with the modelling goal. It should be bear in mind that every additional compartment adds 32 differential-algebraic equations (DAEs) to the model.

Non-Ideal Mixing Kinetics Modelling
The spatial variation of biomass, soluble substrates, pH, specific biogas, and methane production in the CM is discussed below and compared with the CSTR AD of the same volume. Figure 6a shows that the biomass concentration in the CM is 35% less than the corresponding CSTR AD. The reduction in hydraulic retention time (HRT) due to the stagnant zone in the CM affects the growth of biomass; hence, it affects the uptake of soluble substrates and reduces the digester performance. The reduction in the biomass concentration in a CM and the stagnant volume (58% of digester volume) are not linearly proportional because a limited amount of biomass grows in the stagnant zones.
The growth of acetogenic biomass (converting volatile fatty acids to acetate) and methanogenic biomass (converting acetate into biogas) are the two most affected in the CM (Figure 6b). The spatial variation of biomass, soluble substrates, pH, specific biogas, and methane production in the CM is discussed below and compared with the CSTR AD of the same volume. Due to the lower biomass growth in the CM (Figure 6a), the undigested soluble substrates outflow from the digester is higher compared to the CSTR digester. The undigested outflow of the soluble acetate and LCFA is 0.4 g/L and 1.4 g/L, respectively, in the CM compared to only 0.04 g/L and 0.1 g/L in the CSTR AD. Moreover, the specific biogas and methane production are reduced by 14% and 16%, respectively, in the CM compared to the CSTR of the same volume (Figure 6c).
The spatial variation of specific biogas and methane in a CM indicates that all zones of the digester are not equally productive in terms of biogas and methane production. About 88% of specific biogas and methane are produced in zone 2 (C1 and C5), which accounts for only 21% of the digester volume ( Figure 6c). The reason that most of the biogas is produced in zone 2 can be better understood by observing the CFD model velocity contour ( Figure 3) and vectors (Figure 4). The sludge enters zone 2, the outermost compartment ( Figure 5) with a larger diameter and medium velocity; this allows the sludge to spend a longer time in zone 2 and so it is digested before entering the adjacent zones. The outflow sludge from zone 2 (inflow into the remaining zones) contains few undigested substrates. Hence the specific biogas and methane production are very low in the remaining zones (Figure 6b). The specific biogas and methane production is reduced by 15% in the non-ideal mixing (CM) compared to the ideal mixing AD (CSTR AD), which is mainly caused by the ineffective digester volume and limited flow exchange between compartments.  The spatial variation of specific biogas and methane in a CM indicates that all zones of the digester are not equally productive in terms of biogas and methane production. About 88% of specific biogas and methane are produced in zone 2 (C1 and C5), which accounts for only 21% of the digester volume (Figure 6c). The reason that most of the biogas is produced in zone 2 can be better understood by observing the CFD model velocity contour (Figure 3) and vectors (Figure 4). The sludge enters zone 2, the outermost compartment ( Figure 5) with a larger diameter and medium velocity; this allows the sludge to spend a longer time in zone 2 and so it is digested before entering the adjacent As seen in Figure 6d, the biogas composition in zone 2 is 60% CH 4 , and 35% CO 2 , which is similar to several modelled and measured biogas compositions found elsewhere in the literature. In general, the overall biogas composition in the CM and CSTR is not very different. The weighted average of CH 4 and CO 2 composition in the CM is 63% and 31%, respectively, corresponding to 61% and 34% in the CSTR AD.
The pH distribution in the CM ranges between 6.8 to 7.8, while it is 7.0 and constant in the CSTR AD, as indicated in Figure 6e. The pH of zone 2 (C2 and C5), which produces the highest biogas, ranges between 6.8 and 7.2, and this pH range is suitable for the operation of AD. All the remaining zones of the digester pH are greater than 7.2, and it is not a favorable environment for the growth of biomass. The highest pH is 7.8 in the stagnant zone (C8) and is located under the impeller. Figure 6f indicates the concentration distribution of biomass in the CM. The biomass concentration in each compartment is nearly uniform in the AD zones, with a significant velocity magnitude shown in Figure 3. This indicates that the biomass growth is not affected as long as adequate flow exchange between the compartments is maintained. Conversely, in the stagnant zones (C7 and C8), the sludge velocity is nearly zero, and so is the exchange flow rate between the compartments; hence, the biomass concentration is the lowest due to substrate limitation for biomass growth.
In general, the non-ideal kinetic model using a CFD-based CM gives a broad and detailed understanding of the biochemical conversion process of AD. It can be applied for AD and mixer design optimization in the pilot and full-scale AD in which the experimental measurement of biochemical reactions is challenging, and hence the challenges of assuming the digester as a CSTR digester can be detached.

Conclusions
Knowledge-based derivation of the CM and its implementation developed in this work overcomes the assumption of treating the digester in a CSTR. The model gives the detailed local spatial variation of biochemical processes and resolves the limitations of the CSTR digester, having lower computational requirements and taking a shorter time. The findings of this study showed that all digester volumes are not equally important in biogas production as assumed in a CSTR AD, and the performance is reduced by 15% due to the ineffective use of the digester volume. In addition to predicting the non-ideal biogas and methane production, the CM allows the study of the local kinetic behaviors and factors that affect the AD operation, such as pH, which cannot be studied in a simple model such as CSTR. Furthermore, understanding local biokinetic processes using a CFD-based CM leads to the knowledge-based design of AD for better performance, overcoming the effect of non-ideal mixing and reduction in biogas production. Moreover, a CM model removes the challenges of kinetic model stiffness and huge computational demands encountered in integrating the AD kinetics model into the CFD model.

Conflicts of Interest:
The authors declare no conflict of interest.