## Abstract

While computational modelling has increasingly supported wastewater bioreactor engineering, novel numerical techniques have been developed such as the lattice-Boltzmann method (LBM). With vinasse treatment as case study, this work is a first step towards a comprehensive LBM simulator of a continuous-flow anaerobic packed-bed reactor. Extensions from typical models comprise one-dimensional (besides time) dependence, species transport via convection and diffusion, and imposition of either Dirichlet or Danckwerts condition at inlet. The LBM simulator proved to be operational when simulating the bioreactor at different hydraulic retention times (HRTs). Simulated profiles show that stepwise feeding concentrations are smoothed as they are transported towards the bioreactor exit while concentrations increase or decrease in response to generation or degradation kinetics. Good fitting was observed for concentrations of acetic acid (2.1 kg-COD/m^{3} for HRT = 24 h) and butyric acid (1.3 kg-COD/m^{3} for HRT = 16 h) at the exit whereas other concentrations were numerically simulated at proper order of magnitude.

## INTRODUCTION

Computational modelling can provide insights for both fundamental and applied research, particularly when acquiring real-operation or experimental data becomes awkward or even unfeasible (de Souza-Santos 2010). Virtualisation allows one to rapidly examine exploratory scenarios and search for optimal conditions. Accordingly, the shift from empirical (i.e. black-box or data-driven) to physics-based (i.e. white-box or phenomenological) models is vital to benefit from new sensor technology, internet of things and big data (Saguy 2016).

Comprehensive knowledge of bioreactors and wastewater treatment may rely on equations invoking numerical solution methods (Głuszcz *et al.* 2011; Naessens *et al.* 2012). Besides variability in chemical kinetics and composition, mathematical hurdles refer to the mutual interference between fluid velocity and chemical species concentrations (Parco *et al.* 2007). Computational fluid dynamics (CFD) arises as a helpful tool (Brannock *et al.* 2010) and the importance of computational modelling to bioprocess engineering has been recognised while novel numerical methods have been developed and applied (Datta & Sablani 2007).

The lattice-Boltzmann method (LBM) was envisaged in McNamara & Zanetti (1988) as an offspring of another mesoscale method, namely lattice gas cellular automata. LBM has become a powerful technique to numerically simulate bioprocesses and biosystems at different length scales (van der Sman 2007), including bioreactor design (Štumbauer *et al.* 2013). Within the scope of bioreactor simulation, one may apply LBM to computationally model biofilm formation and detachment (Picioreanu *et al.* 2001) as well as biogas bubbles generation and transport (Chen 2010).

From a modelling standpoint, LBM does not directly solve Navier–Stokes equations (Succi 2001) so that it can quite straightforwardly deal with fluid flow, transport phenomena, moving boundaries and multicomponent systems, whether or not coupled to chemical reactions (Sukop & Thorne 2006). On the computational side, LBM is suitable for parallel processing while it renders relatively simpler codes (Mohamad 2011). Those features sound attractive when programming in-house CFD simulators via classic discretisation methods such as finite differences (FDM), finite elements (FEM) or finite volumes (FVM)

LBM avoids the numerical stiffness related to either Galerkin formulations in FEM (Bathe 1996) or pressure-velocity coupling in FVM (Patankar 1980). Discretisation of convective terms (e.g. upwind scheme) may yield false (i.e. numerical) diffusion (Patankar 1980), which is not the case of LBM (van der Sman 2007; Mohamad 2011). False diffusion effects arise in FDM simulations of convective-dominant transport phenomena while they are mitigated by means of LBM (Rabi 2012).

Inspired by the IWA Anaerobic Digestion Model No. 1 (ADM1) (Batstone *et al.* 2002), this work is a first step in the development of a comprehensive in-house LBM simulator of continuous-flow bioreactors for wastewater treatment. While there is no doubt about the efficiency of classic numerical methods and off-the-shelf software, the present work points to LBM as an alternative route to computationally model bioreactors for wastewater treatment.

The physics-based model proposed in this work extends the habitual time-dependent (i.e. dynamic) behaviour to include one-dimensional (1-D) spatial dependence. As far as species transport is concerned, the model considers convection as well as diffusion. The latter has been sometimes ignored as convective-dominant differential equations remain parabolic with respect to spatial dependence (Wu & Hou 2001; Sanchez *et al.* 2004). Consequently, one can dismiss one boundary condition (usually at the exit) while being able to use marching numerical methods such as the Runge–Kutta method (Kreyszig 1993).

Although governing differential equations become elliptic, diffusive transport should be modelled not only when comprehensive knowledge is necessary but also because it is influential in large equipment (Reis-Vasco *et al.* 2000; Gaspar *et al.* 2003). In view of that, convective–diffusive transport has been considered since early versions of our in-house LBM simulators (Durán *et al.* 2015; Okiyama *et al.* 2015; Rabi & Kamimura 2016; Rosa *et al.* 2016). In the present work, the null Neumann condition was imposed at the bioreactor exit while either Dirichlet or Danckwerts boundary conditions were numerically implemented and tested at the bioreactor inlet.

## THEORY

### Physics-based model of continuous-flow bioreactor for vinasse treatment

A cylindrical laboratory-scale anaerobic packed-bed reactor (APBR) for sugarcane vinasse treatment has been taken as case study (Ferraz 2013; Ferraz *et al.* 2014). Vinasse is an effluent from the sugar-ethanol industry whose large-scale use has pointed to ferti-irrigation in sugarcane crops after it has undergone anaerobic treatment (Vlissidis & Zouboulis 1993). Figure 1(a) depicts the experimental APBR assembled within acrylic tubes. Each APBR comprises a feeding module (FM), bed section (BS), effluent collection (EC) area and biogas collection (BC) area, as shown in Figure 1(b).

Low-density polyethylene cylinders (approximate dimensions: length ∼5.0 mm, diameter ∼4.5 mm) were used as supporting medium in the BS module. Those small cylinders randomly filled up the BS module so that bed porosity *ɛ* was observed to range from 0.47 to 0.54 (Ferraz 2013; Ferraz *et al.* 2014). Table 1 summarises length, inner diameter, geometric volume and volume occupied by the liquid phase at each APBR compartment. As the BC module is solely occupied by the gas phase, the volume occupied by the liquid phase could range from *V*_{liq,min} = 2.187 L (for lower porosity) to *V*_{liq,max} = 2.363 L (for higher porosity).

APBR is thermophilic (at 55 °C) and four different volumetric flow rates of vinasse were tested in Ferraz (2013) and Ferraz *et al.* (2014). Flow rates were experimentally set so as to yield hydraulic retention times (HRTs) of 24 h, 16 h, 12 h and 8 h, assessed as . In numerical simulations, corresponding interstitial fluid velocities were evaluated in terms of lower porosity *ɛ*_{min} and bed diameter *d* as . Measured data available for comparisons with numerical simulations (ahead in ‘Results and discussion’ section) refer to species concentrations in the EC outflow, Figure 1(b).

This work proposes a dynamic 1-D model in which coordinate axis *z* is oriented upwards along with vinasse flow. As sketched in Figure 1(c), stratification is assumed so that species concentrations *y _{i}* are functions of time

*t*and coordinate

*z*, i.e.

*y*=

_{i}*y*(

_{i}*z*,

*t*). The 1-D solution domain spans from

*z*= 0 (inlet) to

*z*=

*L*(exit), where length

*L*is numerically set to yield the same HRT from experiments, i.e.

*L*=

*v*

*×*

*HRT*. Table 2 shows the values of APBR operational parameters

*HRT*, ,

*v*and

*L*used in LBM simulations.

*y*are modelled by partial differential equations of the form: where

_{i}*v*is interstitial fluid velocity (which is the same for all species),

*D*is species diffusivity in the fluid phase, and rate includes either source or sink terms concerning species generation and/or consumption through chemical reactions.

_{i}Comprehensive versions of the LBM simulator intend to invoke as many chemical species as those listed in Table 3. For the sake of successful implementation of a preliminary version, this work numerically solves Equation (1) for selected species, namely *y*_{1} to *y*_{5} in Table 3 together with single (say ‘unified’) biomass species *y*_{7} (i.e. rather than invoking species *y*_{71} to *y*_{74}).

### Underlying concepts of LBM

Inspired by the kinetic gas theory, LBM considers any medium (whether solid or fluid) as comprised of fictitious constituent particles following a sequence of synchronised streaming and collision steps in a discrete space (namely a fictitious lattice). During streaming, particles displace from one site to another through links defined by the fictitious lattice assigned to the true medium. As particles simultaneously arrive at lattice sites, they mutually collide so that their discrete velocities are rearranged for subsequent streaming–collision steps. By imposing conservation principles (e.g. mass and momentum) to such repetitive particle dynamics, one may numerically simulate macroscopic medium behaviour (Succi 2001).

The basic mathematical entity in LBM is the particle distribution function giving, at time *t* and about position , the number of particles per unit volume with velocities between and . By taking suitable moments of function *f*, one is able to retrieve macroscopic (i.e. observable) properties such as chemical species concentration, bulk flow velocity, and temperature (Succi 2001; van der Sman 2007; Mohamad 2011).

*f*is governed by Boltzmann's transport equation which, in the absence of external forces, can be written as: where the collision operator Ω gives the variation rate of function

*f*due to collisions between particles. The Bhatnagar–Gross–Krook (BGK) approach was invoked in Equation (2) in order to linearise the collision operator Ω, which means that equilibrium values

*f*

^{eq}of the distribution function are achieved at a rate set by relaxation time Δ

*t*

_{relax}(Qian

*et al.*1992).

LBM is implemented to numerically solve Equation (2) as written in terms of the fictitious lattice, when it becomes known as the lattice-Boltzmann equation (LBE). Particle distribution functions *f _{k}* are then assigned to each streaming velocity

*c*in the lattice. LBM lattices are identified as

_{k}*DnQm*, where

*n*is the problem dimension (e.g.

*n*= 1 for 1-D problem) while

*m*refers to the number of functions

*f*to be solved for each macroscopic property in the true medium.

_{k}The present work proposes a dynamic 1-D model and Figure 2 sketches the basic (and repetitive) linear structure of (a) D1Q2 and (b) D1Q3 lattices. It comprises a central lattice site linked to two neighbouring sites, one in forward (*k* = 1) and the other in backward (*k* = 2) direction, thus leading to streaming velocities *c*_{1} = +*c* and *c*_{2} = −*c*. The difference between D1Q2 and D1Q3 relies on the fact that the latter lattice allows some particles to remain at rest (i.e. null velocity *c*_{0} = 0) after collisions. Usual 2-D and 3-D lattices in LBM are sketched elsewhere (Succi 2001; van der Sman 2007; Mohamad 2011).

*f*(related to each streaming link

_{k}*k*) at axial position

*z*and time

*t*, namely: If Δ

*t*is the advancing time step then streaming velocities are

*c*

_{1}= +Δ

*z*/Δ

*t*and

*c*

_{2}= −Δ

*z*/Δ

*t*. In D1Q3 Equation (3) is written for

*k*= 0, 1 and 2 while in D1Q2 it is solely written for

*k*= 1 and 2, i.e. particle distribution function

*f*

_{0}is disregarded in LBM simulations using D1Q2 lattice.

## NUMERICAL METHOD: LBM SIMULATION OF THE CONTINUOUS-FLOW APBR

As far as LBM simulation of species transport in multicomponent systems is concerned, one may follow two approaches. Species can be simulated either (i) as an active component whose local concentration increase leads to local decrease in background fluid concentration or (ii) as a passive solute carried by the flowing solvent (Sukop & Thorne 2006). In view of LBM simulators already implemented in our research, the present work follows the second approach.

*f*

_{i}_{,k}are assigned to each concentration

*y*invoked in the model. Functions

_{i}*f*

_{i}_{,k}(

*z*,

*t*) become known at any time

*t*and position

*z*by numerically solving Equation (3) written for each species, whose concentration

*y*can then be retrieved as: where the summation runs from

_{i}*k*= 0 to 2 in D1Q3 lattice or from

*k*= 1 to 2 in D1Q2 lattice.

*f*

_{i}_{,k}for all lattice links

*k*are updated from instant

*t*to instant

*t*+ Δ

*t*at all lattice sites (i.e. positions

*z*) in the solution domain. Eventual source or sink terms (e.g. species generation or consumption through chemical reactions) are included at this LBM step so that the following algebraic expression holds: where

*ω*= Δ

_{i}*t*

_{relax,i}/Δ

*t*is known as the relaxation parameter and

*w*are weighting factors fulfilling the condition . Weighting factors for D1Q2 lattice are

_{k}*w*

_{1}=

*w*

_{2}= 1/2 whereas those for D1Q3 lattice are

*w*

_{0}= 4/6 and

*w*

_{1}=

*w*

_{2}= 1/6.

As D1Q2 lattice is used, streaming was implemented in forward (Δ*z*_{1} = +Δ*z*) and backward (Δ*z*_{2} = −Δ*z*) directions. Streaming can be suppressed (with no loss of functionality) for any modelled property explicitly lacking derivatives with respect to spatial coordinates in its governing equation (Durán *et al.* 2015; Okiyama *et al.* 2015; Rabi & Kamimura 2016; Rosa *et al.* 2016).

*t*

_{final}is reached. Besides links set in Equation (4), macroscopic medium is connected to LBM via relaxation parameters

*ω*and equilibrium distribution functions . The former are linked to the transport coefficients so that in the present dynamic 1-D model relaxation parameters

_{i}*ω*are related to species diffusivities

_{i}*D*according to: where Ma and Pe

_{i}*are respectively lattice-based Mach and species-transfer Péclet numbers. Equilibrium distribution functions are linked to bulk fluid (i.e. effluent) velocity according to: where the sign before the lattice-based Mach number is positive for forward (i.e. upward,*

_{i}*k*= 1) streaming or negative for backward (i.e. downward,

*k*= 2) streaming.

*f*

_{i}_{,1}(0,

*t*) is unknown. Danckwerts inlet condition is imposed by combining first-order finite-differences estimate of ∂

*y*/∂

_{i}*z*together with Equation (4), thus leading to: where

*y*

_{i}_{,in}is the species concentration in the feeding effluent. Dirichlet inlet condition can be obtained by letting Pe

*→∞ in the above expression, namely:*

_{i}*f*

_{i}_{,2}(

*L*,

*t*) is unknown. Null Neumann condition is imposed by again combining first-order finite-differences estimate of ∂

*y*/∂

_{i}*z*with Equation (4) so that the final result is: Last but not least, initial conditions for particle distribution functions are implemented as: where

*y*(

_{i}*z*,0) are the initial conditions imposed to chemical species concentrations.

## RESULTS AND DISCUSSION

*y*

_{1}to

*y*

_{5}(in Table 3) together with biomass concentration

*y*

_{7}. In view of that, Table 4 shows corresponding generation/consumption rates as expressed in terms of a mutual auxiliary parameter

*A*, namely:

For the chemical species preliminary invoked, Table 4 also shows initial conditions *y _{i}*(

*z*,0) for Equation (12) as well as feeding values

*y*

_{i}_{,in}for either Danckwerts or Dirichlet inlet conditions, Equation (9) or (10). The latter were measured in Ferraz (2013) and Ferraz

*et al.*(2014) and they refer to typical features of sugarcane vinasse. Null Neumann boundary condition, Equation (11), was imposed to all species at exit.

Besides the operational parameters from Table 2, LBM simulations were performed by setting *k*_{d1} = 0.02 day^{−1}, *Y*_{su} = 0.1, *f*_{C2,su} = 0.41, *f*_{C3,su} = 0.27, *f*_{C4,su} = 0.13, *f*_{H2,su} = 0.19, and *D _{i}* = 0.02 m

^{2}/day. The latter value was a small diffusivity arbitrarily set for all chemical species and a sensitivity analysis is left for future work. LBM parameters were set as Δ

*z*= 0.01 m and Δ

*t*= 0.0005 day in order to render small lattice-based Mach number (namely Ma = 0.1) for numerical stability purposes (Mohamad 2011).

Available experimental data refer to species concentrations *y*_{1} to *y*_{4} measured at EC exit, Figure 1(b), after 1 day of operation (Ferraz *et al.* 2014). For each operation (i.e. HRT) scenario in Table 2, species concentrations numerically simulated at *z* = *L* and *t* = 1 day are compared with aforesaid experimental data in Table 5. Good fitting was observed for some species concentrations (e.g. acetic acid for HRT = 24 h and butyric acid for HRT = 16 h) and other concentrations were numerically simulated at appropriate order of magnitude.

Differences between experimental values and LBM simulations might refer to the fact that chemical parameters have not been fine-tuned. Also, the model must be extended to include further chemical species and corresponding governing equations (together with additional parameters). Those computational modelling tasks are among future developments of the LBM simulator and they fall beyond the scope of this preliminary work.

Species concentrations in Table 5 are in fact a small sample of the numerical results one is able to obtain. The LBM simulator can indeed provide species concentrations *y _{i}* at any axial position

*z*and time

*t*. While those numerical values can be contrasted with measured data, comparisons in this work are limited to experimental values available in Table 5. Yet, one may grasp bioreactor behaviour from simulations and, in what follows, concentration profiles are shown for LBM simulations performed for APBR operated at HRT = 12 h.

For each chemical species invoked in the current preliminary model, Figure 3 shows time-dependent concentrations *y _{i}*(

*L*,

*t*) numerically simulated at the exit. Concentrations simulated under Danckwerts condition at the inlet, Equation (9), are compared to Dirichlet counterparts, Equation (10). Only minor differences are observed when one inlet boundary condition is replaced by the other. In effect, differences are prone to become evident for relatively high diffusivities

*D*(Rabi & Kamimura 2016). LBM simulations suggest that a time lag of approximately 8 h (∼1/3 day) is necessary for concentration changes to become traceable at the bioreactor exit.

_{i}While Figure 3 shows how species concentrations may vary in time at a specific position, the LBM simulator can also depict how concentrations may either spread or cluster within the bioreactor during operation. Accordingly, Figure 4 shows species concentration profiles *y _{i}*(

*z*,

*t*) numerically simulated over the solution domain (i.e. for 0 ≤

*z*≤

*L*) at the following instants:

*t*= 0.25 day (=

*t*

_{final}/4),

*t*= 0.5 day (=

*t*

_{final}/2),

*t*= 0.75 day (=3

*t*

_{final}/4), and

*t*= 1.0 day (=

*t*

_{final}).

From COD profiles simulated at *t* = 0.25 day and *t* = 0.50 day in Figure 4(a), one may note that the stepwise feeding concentration is smoothed while it is transported towards the bioreactor exit. As a result, there is a concentration increase at this position, which is evidenced in Figure 3(a) for *t* < 0.6 day. Subsequently, COD concentration at *z* = *L* is reduced due to its degradation throughout the bioreactor.

As far as acid concentrations are concerned, Figure 4(b)–4(d) suggest that stepwise feeding fronts are equally smoothed as they are transported to the exit. Conversely, those species concentrations are augmented due to acidogenesis over the bioreactor. Similar increasing behaviour was numerically simulated for dissolved H_{2} and biomass as Figure 4(e) and 4(f) respectively suggest, with degradation effects observed in the latter species.

All species concentration profiles in Figure 4 were numerically simulated by imposing null Neumann condition at the exit and Danckwerts condition at the inlet. LBM simulations under Dirichlet inlet condition yielded similar concentration profiles over the solution domain and they are not shown here for brevity. Yet, in Figure 4(f) one may spot the effect of imposing Danckwerts condition as the inlet concentration *y*_{7}(0,*t*) slightly increases over time. This occurs because Danckwerts condition comes from a mass balance accounting for species transport by both convection and diffusion at the inlet (Danckwerts 1958). Conversely, in the case of imposing Dirichlet condition at inlet, concentration *y*_{7}(0,*t*) remains fixed at the feeding value, i.e. *y*_{7,in} = 0.1 kg-COD/m^{3} (as given in Table 4).

Short-term developments of the LBM simulator will refer to (i) including further species in Table 3 (together with governing equations) and (ii) coupling optimisation routines to best-fit model parameters (e.g. fine-tuning diffusivities *D _{i}*) against available experimental data. Long-term improvements will include LBM simulation of fluid (i.e. vinasse) flow while extending the solution domain to two (and eventually three) dimensions.

## CONCLUSIONS

Advances in dynamic modelling of bioreactors were successfully implemented in the LBM simulator. Specifically, this work extends bioreactor simulation to include 1-D spatial (besides time) dependence, species transport through convection as well as diffusion, and imposition of either Dirichlet or Danckwerts condition at the inlet. At its preliminary development stage, the LBM simulator proved to be operational when simulating the continuous-flow APBR at different HRT.

There was good fitting of concentrations of acetic acid (2.1 kg-COD/m^{3} for HRT = 24 h) and butyric acid (1.3 kg-COD/m^{3} for HRT = 16 h) while other species concentrations were numerically simulated at proper order of magnitude. LBM simulations were able to depict in detail how species concentrations change throughout the bioreactor length at specific time instants as well as how they vary in time at the bioreactor exit. Those numerical results provide insights to how chemical species may either spread or cluster within the bioreactor.

## REFERENCES

*.*

*.*

*,*2nd edn.

*PhD Thesis*

*.*7th edn.

*.*

*.*

*.*

*.*

*.*

Anaerobic Digestion of Sugar Cane Vinasse in Acidogenic Fixed Bed Reactor Followed by Methanogenic Reactor Sludge Blanket Type