Characterising the two-phase flow and mixing performance in a gas-mixed anaerobic digester Importance for scaled-up applications

This study aimed to characterise the gas-liquid ﬂ ow and mixing behaviour in a gas-mixed anaerobic digester by improving phase interaction modelling using Computational Fluid Dynamics (CFD). A 2D axisymmetric model validated with experimental data was set up using an Eulerian-Eulerian method. Uncertainty factors, including bubble size, phase interaction forces and liquid rheology were found to signi ﬁ cantly in ﬂ uence the ﬂ ow ﬁ eld. A more reliable and complete validation was obtained by critical comparison and assessment of the referred experimental data, compared to the models reported in other studies. Additionally, justi ﬁ able corrections and predictions in detail were obtained. Mixing was evaluated by trajectory tracking of a large number of particles based on an Euler-Lagrange method. The mixing performance approximated to a laminar-ﬂ ow reactor (LFR) that distinctly deviated from expected continuous stirred tank reactor (CSTR) design, indicating limited enhancement from the applied gas-sparging strategy in the studied digester. The study shows the importance of a proper phase-interaction description for a reliable hydrodynamic characterisation and mixing evaluation in gas-mixed digesters. Validations, bend to experimental data without a critical assessment, may lead to an inaccurate model for further scaled-up applications


Introduction
Stabilisation of biomass and its biochemical conversion into energy-rich biogas are commonly performed in anaerobic digesters, in which good mass transfer and heat transfer rely on a proper mixing.In digesters designed as a continuous stirred tank reactor (CSTR), mixing is commonly achieved by biogas recirculation (Lindmark et al., 2014).However, the mixing is often insufficient in full-scale gas-mixed digesters, due to unaccounted scale-up effects (Bello-Mendoza and Sharratt, 1998;Capela et al., 2009;Terashima et al., 2009).Prevailing short-circuiting or dead zones in practice result in treatment performances below the theoretical potential (Capela et al., 2009;Samstag et al., 2016).In order to optimise mixing, the flow behaviour and mixing performance should be well characterised, and the actual contribution to mixing from the biogas recirculation should be evaluated.Any additional mixing brought forward by the produced biogas resulting from sludge digestion is not taken into account in the design process and the exact impact of evolving biogas has not yet been evaluated.
For research purposes it is very challenging to obtain a clear insight of the mixing behaviour in an opaque and gas tight system when using only experimental approaches.A general mixing pattern can be determined using tracer tests (Capela et al., 2009;Smith et al., 1993;Terashima et al., 2009), but local flow fields responsible for any poor mixing cannot be determined from the results of tracer residence time distribution (RTD).Advanced noninvasive techniques, such as computer automated radioactive particle tracking (CARPT) and computed tomography (CT) can provide detailed information of flow fields (Karim et al., 2004;Varma and Al-Dahhan, 2007), but are yet not applicable to full scale facilities.
In the gas-mixed digesters, sludge flow and mixing are mainly driven by energy from the gas-sparging.Therefore, a proper characterisation of the phase interaction between gas and liquid (sludge) is crucial.However, in the aforementioned studies, the two-phase flow modelling was still limited with considerable simplifications and uncertainties.Some studies eliminated the twophase interaction region and essentially simulated a single-phase flow (Coughtrie et al., 2013;Karim et al., 2007).Although sludge rheology has been recognised as important for design and operation of sludge treatment processes (Eshtiaghi et al., 2013), no rheological measurements of the used sludge were reported (Karim et al., 2007).Hence, different types of liquids were used, including Newtonian (Coughtrie et al., 2013;Vesvikar and Al-Dahhan, 2005) and non-Newtonian fluids (Karim et al., 2007;Wu, 2010a), to validate the same experimental data (Karim et al., 2004).Despite the different used rheology, similar good agreements were reported.Moreover, the model validations were incomplete involving only a part of the experimental results, thus not convincing enough.Although a more detailed model validation was done by Dapelo et al. (2015), the studied reactor was essentially a bubble column with a bottom-mounted nozzle (Dapelo et al., 2015;Dapelo and Bridgeman, 2018), different from the full scale digesters with vertical-hanging gas lances.
Therefore, the aforementioned description of the two-phase flow and limited model validation, should be improved.This study aimed at developing a reliable CFD model in a scaled-down digester, in order to further guide process assessment and optimisation in scaled-up systems.The focus was on modelling the gassludge interaction process, and detailed model validation and mixing evaluations were carried out.

Computational domain and mesh
A computational domain was created based on the gas-lift digester studied by Karim et al. (2004), which was a 7.2 L cylindrical tank concentrically mounted with a gas injection pipe (0.5 cm inner diameter) and a draft tube (4.4 cm inner diameter) (Karim et al., 2004).More detailed information is given in Fig. 1A.Considering the geometry and layout, a 2D axisymmetric domain was developed.As shown in Fig. 1B, a central axis is used, a free surface for the liquid phase is defined at the top, and the other boundaries are constructed according to the original design.
Structured grids were used for meshing, and finer grids were applied within the sensitive regions of phase interaction around the gas pipe and the draft tube (Fig. 1B).To check grid independency, 5 characteristic grid sizes (A through E) from coarse to fine were created, which are shown in Table 1.

Governing equations and models
Considering the focus on bubbly flow and not on single bubble motion, an Eulerian-Eulerian model was utilised, in which the gas phase is treated as a continuum, similar to the liquid phase.For the sake of simplicity, solid particle motion, interfacial mass transfer and heat transfer were neglected, and a single bubble diameter set for the gas phase.Transient simulations were implemented.Hence, the governing equations given below were solved for each phase.

Volume fraction equation
The volume of each phase in a control cell is determined as (1) where V denotes the volume, a the volume fraction, n the number of phases, and subscript i refers to the phase (liquid or gas).

Mass conservation equation
v vt where r denotes the density, and v . denotes the velocity.

Momentum conservation equation
v vt where p denotes the pressure, t the stress-strain tensor, and g . the gravitational acceleration.F . PI denotes the phase interaction forces between the gas and the liquid per unit volume, in which drag force F .D , lift force F .L and turbulent dispersion force F .T are taken into account.Specific models were applied to close these phase interaction forces.

1) Drag force
The drag force is represented by solving the interphase momentum exchange term from the gas phase g to the liquid phase l, as expressed in Eq. (6).
where t g denotes the particulate relaxation time, d g the diameter of the secondary (gas) phase, m l the dynamic viscosity of the primary (liquid) phase, and A in the interfacial area per unit mixture volume.
The drag function f D is solved by the Schiller and Naumann Model (Naumann and Schiller, 1935).
The lift force stems from velocity gradients of the primary phase, and is expressed as where C L denotes the lift coefficient and is defined by the Legendre-Magnaudet model (Legendre and Magnaudet, 1998).
where C L, lowRe , C L, highRe and Sr denotes the parameters in the model: and Re u the vorticity Reynolds number: 3) Turbulent dispersion force The turbulent dispersion force is related to the effect of interphase turbulent momentum transfer, solved using the Lopez de Bertodana Model (de Bertodano, 1991).
where C T denotes a constant (generally set to 1), and k l the turbulent kinetic energy of the liquid phase, which is calculated by solving the k-ε model.

Turbulence closure
Regarding the reported considerable turbulent kinetic energy (Karim et al., 2004), turbulence needed to be modelled.In this process, turbulence originates from the large velocity differences induced by gas-sparging, which is strongly correlated to bubbly regions.Moreover, modified RANS models have been found to be appropriate for a low-Re non-Newtonian flow (Wu, 2010b).Considering the localised bubbly and turbulent distribution, the  (Karim et al., 2004), and (B) the computational domain with boundaries and typical mesh in this study.

Table 1
Mesh setting by 5 characteristic grid sizes.RNG k-ε Mixture Turbulence model was applied, which is more suitable to (nearly) stratified and low-Re multiphase flows (Orszag et al., 1993).
where a k and a ε are the inverse effective Prandtl numbers for k and ε, respectively (a k ¼ a ε z1.393 for high Re).In this model, mixture (denoted by m) properties of both phases are involved and are computed from Eq. ( 18).
where G k,m denotes the production of k.The Differential Viscosity Model was applied to better solve low-Re flows, so the effective The parameter C 2ε * is calculated by Eq. ( 20).
The other constants are shown below: Due to the finer grids, the first grid layer attaching wall boundaries could be in the viscous sub-layer of turbulence (Y þ <5).Hence, an Enhanced Wall Function was applied to correctly calculate the near-wall region flow field covering both linear and logarithmic patterns.

Phases
The detailed phase physical properties are shown in Table 2.For the primary phase, different fluids were used, including Newtonian (water) and non-Newtonian (sludge 1 and 2) that were described by the power law.
where K denotes the flow consistency index, _ g the shear rate, and n the flow behaviour (power-law) index.Sludge 1 data were from literature (Achkari-Begdouri and Goodrich, 1992).Sludge 2 was the waste activated sludge sampled from Wastewater Treatment Plant De Groote Lucht (The Netherlands); the rheology was measured according to our previous work (Wei et al., 2018), and total solids (TS) were measured according to standard methods (APHA, 2012).The secondary phase was set as spherical air bubbles, with a diameter in the range of 1e10 mm.

Euler-Lagrange method
An Euler-Lagrange method was applied to track particle trajectories and evaluate the mixing.Spherical particles with the same settings to the experiment (Karim et al., 2004) were used in a oneway coupling simulation.The turbulence dispersion (by Discrete Random Walk model), virtual mass and pressure gradient effects on particles were considered.Thus the particle motion is solved by: where the subscript pt denotes the particles, v l the mean velocity, z a normally distributed random number, and t pt the particle relaxation time.The virtual mass force (F VM ), pressure gradient force (F PG ) and Saffman lift force (F SL ) in unit volume are solved by: where C VM is the virtual mass coefficient (¼ 0.5), and n l the liquid kinematic viscosity.A larger number of particles (9063) were released and tracked at different gas flow rates.

Boundary conditions
Involved boundary conditions included axisymmetric centreline, velocity inlet, degassing top and wall.The axisymmetric centreline was set for creating a 3D axisymmetric circumstance.The velocity inlet was for gas injection.The degassing top was for discharging the gas phase while retaining the liquid phase.The wall condition including no-slip for liquid and free-slip for gas was set at the other boundaries.

Main simulation settings
In the transient simulations, time step setting obeyed the Courant-Friedrichs-Lewy (CFL) condition and the First Order Implicit scheme was used.The Phase Coupled SIMPLE scheme was utilised for pressure-velocity coupling.High order schemes including Second Order Upwind, QUICK (Leonard and Mokhtari, 1990) and MUSCL (van Leer, 1979) were set for improving simulation accuracy.Regarding the relatively complicated process modelled and performance achieved, the simulations were done in a stringent way with satisfactory convergence: 1) a desired criterion of 1 Â 10 À5 was used for residuals of most parameters, while a more reasonable order of 10 À4 was used for the continuity and gas velocity residuals; because residual values of the aforementioned parameters were difficult to reach 1 Â 10 À5 , even using a quite small time step of 1 Â 10 À6 s (1 Â 10 À4 s already fits the CFL condition).2) the gas mass flow rate and liquid flux were also monitored and achieved satisfactory balance (imbalance < 5%).The simulations were implemented using the commercial package ANSYS-Fluent 16.2 mainly on a Dell Optiplex 7010 computer, with Intel Core i5-3740 and 8 GB RAM.When running two parallel processors, the computational time ranged between 1 and 10 days.

Grid independence test
A grid independence test was performed first, using the Grid Convergence Index (GCI) method (Celik et al., 2008;Roache, 1994).
The axial velocities at heights of 11.25 cm and 18.25 cm were used for the assessment.Mesh D and E were found to have quite similar profiles and values of the velocity; and the GCI results (GCI CD : 3.3 ± 3.5%, GCI DE : 0.7 ± 1.1%) were in an asymptotic range close to 1 (1.02 ± 0.02), indicating that the grid independence was achieved.Considering smaller number of grid cells and less computational consumption, Mesh D was selected as optimal for the following simulations.The detailed results and interpretation can be found in the Supplementary Material (Fig. S1 and Table S1).

Bubble size
As mentioned in the Introduction, some uncertainties needed to be investigated.Although no experimental data are available to support, proper identification of a characteristic bubble size is important in the model calculations.To achieve a clear estimation, only the drag force was applied for phase interaction, and two distinct fluids water (Newtonian) and sludge 1 (non-Newtonian (Achkari-Begdouri and Goodrich, 1992)) were used.
According to the references (Lehr et al., 2002;Polli et al., 2002;Sch€ afer et al., 2002;Vesvikar and Al-Dahhan, 2005;Wu, 2010a), bubbles with diameters ranging from 1 mm to 10 mm can be expected.Fig. 2 shows the liquid axial velocity profiles at 1/2 depth, zooming in inside the draft tube.Although the peak became smaller, similar velocity profiles of the water were obtained as the bubble size increased, especially larger than 3 mm (Fig. 2A).However, different profiles were obtained in sludge 1.As the bubble size increased, a more substantial decrease was found in the peak value, and the peak position also moved further from the centreline (Fig. 2B).Similar tendencies were found at the other heights in both fluids (data not shown).
The velocity independency on the bubbles over 3 mm agreed with a previous study using water as well (Vesvikar and Al-Dahhan, 2005), in which the corresponding bubble size was 2 mm, and 10 mm was set in their simulations.However, the velocity was quite sensitive to the bubble size in sludge 1.Hence, the bubble size setting (1 mm) in the study of Wu (2010a) using the same sludge data needs to be evaluated; in experimental measurements, the characteristic bubble sizes in free air-water bubbly flow were reported from 1 to 5 mm (Lehr et al., 2002;Polli et al., 2002;Sch€ afer et al., 2002).Moreover, the bubble size was found to increase with the liquid apparent viscosity (Sch€ afer et al., 2002).Therefore, the bubble size settings in previously reported simulations were not suitable: 10 mm (Vesvikar and Al-Dahhan, 2005) seemed too big for water and 1 mm (Wu, 2010a) too small for sludge.Considering the typical bubble sizes in reality and the apparent viscosity influence, 5 mm was set in the following simulations.

Phase interaction forces
Presence of the pipe and the draft tube in the studied digester created more changes in the local liquid velocity than that in conventional bubble columns; as mentioned in Section 2.2, the bubbly process induces high velocity differences and fluctuations.Hence, some potential influences related to velocity gradient and  turbulence need to be considered.Compared to the previous studies (Coughtrie et al., 2013;Vesvikar and Al-Dahhan, 2005;Wu, 2010a), more phase interaction forces, including lift force and turbulent dispersion force, were involved in this study.
As shown in Fig. 3A and B, with the addition of the lift force and/ or the turbulent dispersion force, both water and sludge 1 get a lower but wider velocity peak inside the draft tube and a stronger flow outside, even though the liquid phase got low Re (10 2 ).The discrepancy was also reflected in gas phase distribution results, in which the gas phase was distributed over a much larger domain when applying all the forces (Fig. 3C).Hence, the results indicated important roles for the two forces, leading to broaden the distribution of the gas phase and the liquid flow.
Unlike buoyancy and the drag force causing vertical rising, the lift force accounted for horizontal bubble motion.As shown in Fig. the sharp velocity gradient in the near-pipe region (r/ R < 0.05) can result in a considerable bubble dispersion and turbulence, which further induces more liquid flows.However, the more viscous sludge 1 weakened this effect, leading to a smaller velocity change than that of the water.Furthermore, a similar tendency was found by applying larger bubbles and only the drag force in sludge 1 (Fig. 2B), which might result from the shear-ratedependent rheological behaviour.Hence, for a non-Newtonian fluid, the high velocity differences not only induced more phase interaction forces, but also led to localised differences in the apparent viscosity, making the two-phase flow more complicated.

Liquid apparent viscosity
Sludge 2 measured in this study was also used to investigate the impact of apparent viscosity.Regarding the reported 4.5% (Karim et al., 2004), sludge 1 and 2 had comparable TS concentrations, but distinct rheological behaviour, due to the sludge origin (Table 2).
As shown in Fig. 4, the predicted velocity is quite different: sludge 1 shows a more constrained flow than the water, demonstrated by a narrowed velocity profile inside the draft tube, and a smaller vortex outside with a separation point (the liquid velocity is zero and changes direction) closer to the centreline.Moreover, an extremely constrained flow was found for the much more viscous sludge 2 without the induced vortex or flow.The flow field change also indicated a trend to a highly concentrated and unrealistic gas phase distribution, implying the limitation of the Eulerian-Eulerian method for extremely viscous fluids.
These results were different from previously reported results, in which it was described that a rheological change did not affect mixing behaviour significantly (Karim et al., 2007).It should be noted that in the latter study two fluids were applied that had comparable apparent viscosities (pseudoplastic: 2e9 mPa s and dilatant: 8e32 mPa s ( _ g range 1e1000 s À1 )).This is not large enough to cover the rheological differences found in practice.As illustrated in our study, a change in sludge rheology could lead to a quite different flow pattern under the same gas-sparging conditions.It can be concluded that accurate characterisations of sludge rheology are crucial for proper flow prediction and mixing evaluation.

Validation and further prediction
Considering the aforementioned incomplete validation in the previous research (Coughtrie et al., 2013;Karim et al., 2007;Vesvikar and Al-Dahhan, 2005;Wu, 2010a), and the sensitive impacts on the flow field found in Section 3.2, the validation of the proposed model was improved.

Velocity field
Compared to the experimental velocity field in the whole domain (Karim et al., 2004), general flow patterns were in agreement, including the vortex, the major recirculation loop, and the strong/weak flow locations.However, consistent with Section 3.2, discrepancies were considerable between the water and sludge 1. Highly similar to the previous study using water as well (Vesvikar and Al-Dahhan, 2005), as shown in Fig. 5, the water has a big vortex covering over one half of the whole depth, and the vortex centre is at 17.5 cm height, lower than the top of the draft tube.However, the measured vortex was smaller, covering one third of the whole depth and having a higher centre (19.25 cm height).Hence, this considerable overestimation of the flow cannot be used to conclude a good agreement as reported (Vesvikar and Al-Dahhan, 2005).Although predicting a slightly constrained flow field with a smaller vortex closer to the centreline, sludge 1 showed a better approximation.It indicated that the unreported apparent viscosity of the experimental sludge might have been between that of the water and sludge 1, even better resembling sludge 1.
The axial velocity profiles of sludge 1 at three heights labelled in Fig. 5 are shown in Fig. 6.Considering the uncertain settings discussed above, the overall agreement was satisfactory.For all gas flow rates, the velocity outside the draft tube (r/R > 0.22) only showed small differences with the almost fixed separation point (0.64) that was smaller than the measured 0.71 (Karim et al., 2004).It confirmed the slightly constrained flow mentioned above, and also indicated the difficulty to further enhance this part of the flow under the given gas-sparging conditions.However, more discrepancies were found inside the draft tube, which are further discussed in Section 3.3.3.

Gas phase distribution
Besides velocity field, gas holdup is another important indicator that has not been validated before.Most bubbles were observed to rise inside the draft tube and were mainly confined in the 1.8 mm radius central region (Karim et al., 2004), while almost no gas was detected outside the draft tube (Karim et al., 2007).As shown in Fig. 7A, the gas phase plume, containing a cluster with high volume fractions (>0.2), is consistent with the descriptions above, and with other observations about separated rising bubbles without coalescence (Dapelo et al., 2015;Vesvikar and Al-Dahhan, 2005).In Fig. 7B, the predicted maximum gas volume fraction of 0.53 at 20 cm height is comparable to the experimental value of 0.48 (Karim et al., 2004), which is much better than the value (<0.02) reported previously (Vesvikar and Al-Dahhan, 2005).
Although a similar tendency with a small gas volume peak was predicted just above the draft tube (r/R ¼ 0.21), the simulation demonstrated a more confined and highly underestimated gas distribution in the other parts of the reactor.To assess this large deviation, the experimental data accuracy had to be evaluated.The measured peak (at 1.8 mm) was positioned inside the pipe (radius 3.2 mm, the dash line in Fig. 7B), which seems incorrect.Additionally, a considerable portion of the gas holdup outside the draft tube was reported to result from data discontinuity (Karim et al.,Fig. 4. Velocity of water, sludge 1 and sludge 2, profiles at heights of (A) 11.25 cm, (B) 18.25 cm, and (C) contours in the whole domain; gas flow rate 28.32 L/h.Fig. 5. Simulated and experimental (Karim et al., 2004) velocity fileds in the digester, gas flow rate 84.96L/h.2004).Hence, the deviation could have also originated from inaccurate data collection.
A conserved gas flux calculation (gas flux per cell ¼ gas volume fraction Â gas velocity) was performed for the simulation data, but couldn't be done for the experimental results due to lack of reported gas velocity data.The aforementioned rare bubble occurrence outside the draft tube correlated to a quite small effective gas velocity in the studied period, leading to a negligible effective contribution to the total gas flux.Hence, the dominant part of this balanced gas flux was still from the central region where most bubbles occurred and where the predicted and experimental data were close to each other.Additionally, the realistic random motion of a single bubble is essentially difficult to be represented by the applied Eulerian-Eulerian method with a RANS model, also resulting in the constrained gas phase distribution.Therefore, the large deviation in Fig. 7B did not mean a failure of the simulation; contrarily, the simulation results of the conserved gas phase had a peak in the key region and a bubbly pattern comparable to the experimental data, still indicating a good agreement.

Critical validation and further prediction of the flow field
Considering that the sludge used in the experiment was expected to be a little less viscous than sludge 1, the velocity deviation inside the draft tube requires some further discussion.
To assess data reliability, the flux balance of the liquid axial velocity was evaluated.As shown in Table 3, the simulation results were excellent with a mean imbalanced flux of 0.83 ± 0.75%; but the experimental data showed a large value of 11.3 ± 9.4%.The only measured dataset (in Fig. 7B) even got higher imbalanced flux values (in parentheses in Table 3), also implying that the measured gas holdup may not be accurate enough.In addition, the velocity data were not consistent in different studies, although from the same measurements (Karim et al., 2004;Vesvikar and Al-Dahhan, 2005).As shown in Fig. 6, the experimental velocity profiles in D showed more discrepancies with A-C, having peaks closer to the centreline and a different profile above the draft tube.Moreover, Fig. 6AeC shows considerable velocity values at the position of the draft tube at 11.25 cm height, which does not obey the theoretical zero-velocity assumption on the no-slip walls.The large standard deviation bars of the velocity at the same height in Fig. 6D demonstrated strong fluctuations.Therefore, the accuracy and consistency of the experimental velocity data were not good enough, which was affected by the strong flow fluctuation and limited resolution inside the draft tube (only 4 data points).
Unlike the previous research (Vesvikar and Al-Dahhan, 2005;Wu, 2010a), the validation in this study was not assessed only by Fig. 6.Axial velocity profiles of sludge 1 and measurements (Karim et al., 2004) at 3 heights, gas flow rates of (A) 28.32 L/h, (B) 56.64 L/h, and (C) 84.96 L/h, and (D) the results of Vesvikar and Al-Dahhan (2005).
the degree of overlap to the experimental data.Even achieving a high agreement, the obtained model may not be optimised if having no well-balanced mass/flux and some uncertain parameter settings.Hence, adjusting the model to obtain a perfect fit to the experimental data could be risky for any further application of the developed model, which must be carried out with care.The aforementioned deviation may also be related to the nature of the Eulerian-Eulerian method for the model.In the bubble-driven process, the bubbles-concentrated region, represented by a high gas volume fraction, can have a strong buoyancy effect.The generated large velocity difference across the gas-liquid interface can then induce strong shear and fluctuation of the liquid velocity.Hence, the positions of the high gas phase fraction (a G ), the high liquid velocity (v L ), and the strong momentum transfer (F PI,L ) are correlated.On one hand, measurements with 4 points may be too limited to characterise the quite localised and unstable bubbleinduced flow inside the draft tube.On the other hand, the assumption of a fluid instead of separate bubbles for the gas phase, makes the Eulerian-Eulerian with RANS model not specific enough to simulate a bubble-induced non-Newtonian flow.However, the improved two-phase simulation presented in this study not only achieved good agreement with the experimental data, but also made reasonable corrections and further predictions of the hydrodynamics in detail.Hence, the elaborated scaled-down model could be considered for further scale-up applications.

Mixing evaluation
Fig. 8A shows the typical particle trajectories in the simulation, and various paths are obtained in the particles labelled by different colours.Most formed big loops back into the draft tube, which agreed with the experimental data (Karim et al., 2004) shown in Fig. 8C.As shown in Fig. 8B, the typical instantaneous particle distribution is not homogeneous, but confined to form some major paths.Some regions, including the big vortex, bottom and corner of the digester, showed almost no particles travelling.Mixing performance and RTD were mainly evaluated using particle recirculation time, defined as the time for a particle to return at the released position.Totally 9063 particles were tracked while applying three gas flow rates, and normalised particle number and time (q) were applied, which referred to the total particle number and theoretical mean recirculation time, respectively.The ideal models of CSTR, plug-flow reactors (PFR) and laminar-flow reactors (LFR) were also used for evaluation.
As shown in Fig. 9A, the major peaks are almost independent from the gas flow rates.However, one separated smaller peak (over q ¼ 1) was found at 28.32 and 56.64 L/h, while a much smoother trend was found at 84.96 L/h.The results of key indicators for RTD and mixing in Table 4 showed similar short-circuiting indexes (q 10 and q p ) smaller than the ideal models (0.5 for LFR and 1.0 for PFR), and a considerable increase in the mixing indexes (q 90 , q g , Mo, Th and s 2 ) when the gas flow rate increased.This indicates the occurrence of short-circuiting and dead zones, part of which could be mitigated by increasing the gas flow rate, as can be seen from the decrease in number of small peaks.However, the enhancement of mixing and dispersion was limited, reflected by the limited change of the major particle recirculation time at increased gas flow rates.
The mixing behaviour was further characterised using cumulative distribution results, combined with the ideal models.As shown in Fig. 9B, besides the steep increase in the beginning, the curves of  (Karim et al., 2004); gas flow rate 28.32 L/h.

Table 3
Imbalanced flux (%) of the axial velocity at the three heights, results in parentheses were based on the measured gas holdup data (Karim et al., 2004) 28.32 and 56.64 L/h are approaching to the ideal PFR curve after the ideal mean recirculation time (q ¼ 1), due to the considerable number of particles unable to recirculate (Table 4: 8.7% and 5.3%, respectively).However, the trapped particles in 84.96 L/h were much fewer (0.7%) and the curve approached the ideal LFR.Therefore, in the given system, the contribution from the gassparging to mixing enhancement was more related to the convection.When the sludge exhibited non-Newtonian and viscous behaviour, the diffusion and dispersion from turbulence became weakened, leading to a more plug-flow pattern with constrained flow regions (short-circuiting), which was demonstrated by the belt-like distributions with condensed particles in Fig. 8B.Increasing the gas flow rate could mitigate the problem, but could not change the intrinsic poor mixing regions, and the enhancement was limited to an LFR pattern that is not designed.Hence, to operate the digester in the desired CSTR mode, the limitations of the gassparging strategy should be concerned, and improvements for the commonly used gas-lance layout are necessary for scaled-up systems.

Conclusions
Bubble size, phase interaction forces, and liquid rheology significantly impacted the two-phase flow hydrodynamics in the studied anaerobic digester.Predictions of the gas phase dispersion and liquid velocity relied on characterisation of the phase interaction forces, which mainly depended on the applied gas-sparging strategy.Accurate model development for scaled-up applications requires a critical assessment of experimental data that are used for model validation.Evaluation of the mixing showed a performance that was approximated with a laminar-flow reactor (LFR), distinctly deviating from the expected CSTR design.Results indicated that the applied gas-sparging strategy for mixing is not very effective for reaching a CSTR flow regime.The results underline the importance of a proper phaseinteraction description for achieving a reliable hydrodynamic characterisation and mixing evaluation in gas-mixed digesters.

Fig. 1 .
Fig. 1. (A) Digester geometry(Karim et al., 2004), and (B) the computational domain with boundaries and typical mesh in this study.

Fig. 2 .
Fig. 2. Liquid axial velocity profiles of different bubble sizes inside the draft tube in (A) water and (B) sludge 1 at 11.25 cm height; gas flow rate 28.32 L/h.

Fig. 3 .
Fig. 3. Results of applying different phase interaction forces: axial velocity profiles of (A) water and (B) sludge 1 at 18.25 cm height, and (C) gas phase distribution in sludge 1; gas flow rate 28.32 L/h.

Fig. 7 .
Fig. 7. (A) Predicted gas phase distribution in the whole domain, and (B) at 20 cm height combined with experimental data(Karim et al., 2004); gas flow rate 28.32 L/h.

Fig. 9 .
Fig. 9. Normalised (A) RTD curve and (B) cumulative distribution curve of the particles in the three gas flow rates.

Table 2
Physical properties of the applied phases. .

Table 4
Indicators of the particles at the three gas flow rates for RTD and mixing evaluation, normalised time applied.