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/m3 for HRT = 24 h) and butyric acid (1.3 kg-COD/m3 for HRT = 16 h) at the exit whereas other concentrations were numerically simulated at proper order of magnitude.

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.

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

Figure 1

Laboratory-scale cylindrical APBR for vinasse treatment: (a) picture of experimental bioreactors; (b) sketch comprising feeding module (FM), bed section (BS), effluent collection module (EC) and biogas collection module (BC), adapted from Ferraz (2013); (c) time-dependent one-dimensional approach for BS, where ε is bed porosity.

Figure 1

Laboratory-scale cylindrical APBR for vinasse treatment: (a) picture of experimental bioreactors; (b) sketch comprising feeding module (FM), bed section (BS), effluent collection module (EC) and biogas collection module (BC), adapted from Ferraz (2013); (c) time-dependent one-dimensional approach for BS, where ε is bed porosity.

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 Vliq,min = 2.187 L (for lower porosity) to Vliq,max = 2.363 L (for higher porosity).

Table 1

Length, diameter, geometric volume and volume occupied by liquid phase of each APBR module

DimensionFeeding moduleBed section
Effluent collectionBiogas collection
Porosity = 0.47Porosity = 0.54
Length (m) 0.100 0.500 0.500 0.100 0.050 
Diameter (m) 0.080 0.080 0.080 0.080 0.080 
Geometric volume (10−3 m3 = L) 0.503 2.513 2.513 0.503 0.251 
Volume with liquid (10−3 m3 = L) 0.503 1.181 1.357 0.503 – 
DimensionFeeding moduleBed section
Effluent collectionBiogas collection
Porosity = 0.47Porosity = 0.54
Length (m) 0.100 0.500 0.500 0.100 0.050 
Diameter (m) 0.080 0.080 0.080 0.080 0.080 
Geometric volume (10−3 m3 = L) 0.503 2.513 2.513 0.503 0.251 
Volume with liquid (10−3 m3 = L) 0.503 1.181 1.357 0.503 – 

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 yi are functions of time t and coordinate z, i.e. yi = yi(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.

Table 2

Values of APBR operational parameters to perform LBM simulations

HRT (day) (m3/day)v (m/day)L (m)
0.002362 1.000 1.000 
2/3 0.003544 1.500 1.000 
1/2 0.004725 2.000 1.000 
1/3 0.007087 3.000 1.000 
HRT (day) (m3/day)v (m/day)L (m)
0.002362 1.000 1.000 
2/3 0.003544 1.500 1.000 
1/2 0.004725 2.000 1.000 
1/3 0.007087 3.000 1.000 
In view of model comprehensiveness and scale-up towards large equipment, convective–diffusive transport is considered in the fluid phase. Hence, chemical species concentrations yi are modelled by partial differential equations of the form:
formula
(1)
where v is interstitial fluid velocity (which is the same for all species), Di is species diffusivity in the fluid phase, and rate includes either source or sink terms concerning species generation and/or consumption through chemical reactions.

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 y1 to y5 in Table 3 together with single (say ‘unified’) biomass species y7 (i.e. rather than invoking species y71 to y74).

Table 3

Species to be considered in comprehensive versions of the LBM simulator of continuous-flow APBR

yiConcentration ofUnitsyiConcentration ofUnits
y1 Chemical oxygen demand (COD) kg-COD/m3 y8 Acetate (ethanoate) ion M = kmol/m3 
y2 Acetic (ethanoic) acid kg-COD/m3 y9 Propionate (propanoate) ion M = kmol/m3 
y3 Propionic (propanoic) acid kg-COD/m3 y10 Butyrate (butanoate) ion M = kmol/m3 
y4 Butyric (butanoic) acid kg-COD/m3 y11 Hydrogen ion (H+M = kmol/m3 
y5 Dissolved hydrogen gas (H2kg-COD/m3 y12 Inorganic carbon kmol-C/m3 
y6 Ethanol (ethyl alcohol) kg-COD/m3 y13 Dissolved carbon dioxide kmol-C/m3 
y71 Biomass for sugar degradation kg-COD/m3 y14 Bicarbonate ion kmol-C/m3 
y72 Biomass for H2 degradation kg-COD/m3 y15 Carbon dioxide in the biogas kmol-C/m3 
y73 Biomass for propionic acid degradation kg-COD/m3 y16 Hydrogen gas (H2) in the biogas M = kmol/m3 
y74 Biomass for butyric acid degradation kg-COD/m3    
yiConcentration ofUnitsyiConcentration ofUnits
y1 Chemical oxygen demand (COD) kg-COD/m3 y8 Acetate (ethanoate) ion M = kmol/m3 
y2 Acetic (ethanoic) acid kg-COD/m3 y9 Propionate (propanoate) ion M = kmol/m3 
y3 Propionic (propanoic) acid kg-COD/m3 y10 Butyrate (butanoate) ion M = kmol/m3 
y4 Butyric (butanoic) acid kg-COD/m3 y11 Hydrogen ion (H+M = kmol/m3 
y5 Dissolved hydrogen gas (H2kg-COD/m3 y12 Inorganic carbon kmol-C/m3 
y6 Ethanol (ethyl alcohol) kg-COD/m3 y13 Dissolved carbon dioxide kmol-C/m3 
y71 Biomass for sugar degradation kg-COD/m3 y14 Bicarbonate ion kmol-C/m3 
y72 Biomass for H2 degradation kg-COD/m3 y15 Carbon dioxide in the biogas kmol-C/m3 
y73 Biomass for propionic acid degradation kg-COD/m3 y16 Hydrogen gas (H2) in the biogas M = kmol/m3 
y74 Biomass for butyric acid degradation kg-COD/m3    

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

Particle distribution function f is governed by Boltzmann's transport equation which, in the absence of external forces, can be written as:
formula
(2)
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 feq of the distribution function are achieved at a rate set by relaxation time Δtrelax (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 fk are then assigned to each streaming velocity ck in the lattice. LBM lattices are identified as DnQm, where n is the problem dimension (e.g. n = 1 for 1-D problem) while m refers to the number of functions fk to be solved for each macroscopic property in the true medium.

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 c1 = +c and c2 = −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 c0 = 0) after collisions. Usual 2-D and 3-D lattices in LBM are sketched elsewhere (Succi 2001; van der Sman 2007; Mohamad 2011).

Figure 2

Basic (and repetitive) linear structure of (a) D1Q2 and (b) D1Q3 lattices to perform LBM simulations in one-dimensional solution domains. Null velocity c0 = 0 (after collisions) is solely allowed in D1Q3 lattice.

Figure 2

Basic (and repetitive) linear structure of (a) D1Q2 and (b) D1Q3 lattices to perform LBM simulations in one-dimensional solution domains. Null velocity c0 = 0 (after collisions) is solely allowed in D1Q3 lattice.

LBE-BGK is written for each particle distribution function fk (related to each streaming link k) at axial position z and time t, namely:
formula
(3)
If Δt is the advancing time step then streaming velocities are c1 = +Δzt and c2 = −Δzt. 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 f0 is disregarded in LBM simulations using D1Q2 lattice.

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.

Accordingly, particle distribution functions fi,k are assigned to each concentration yi invoked in the model. Functions fi,k(z,t) become known at any time t and position z by numerically solving Equation (3) written for each species, whose concentration yi can then be retrieved as:
formula
(4)
where the summation runs from k = 0 to 2 in D1Q3 lattice or from k = 1 to 2 in D1Q2 lattice.
Space–time discretisation of Equation (3) renders an algebraic equation whose computational evolution is implemented in two iterative steps. Time evolution is accomplished during the collision step so that particle distribution functions fi,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:
formula
(5)
where ωi = Δtrelax,it is known as the relaxation parameter and wk are weighting factors fulfilling the condition . Weighting factors for D1Q2 lattice are w1 = w2 = 1/2 whereas those for D1Q3 lattice are w0 = 4/6 and w1 = w2 = 1/6.
Space evolution is achieved in the streaming step as collision outcomes are transmitted to neighbouring lattice sites in all streaming directions as:
formula
(6)

As D1Q2 lattice is used, streaming was implemented in forward (Δz1 = +Δz) and backward (Δz2 = −Δ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).

LBM simulation consists of iterating Equations (5) and (6) throughout the solution domain until the final simulation instant tfinal is reached. Besides links set in Equation (4), macroscopic medium is connected to LBM via relaxation parameters ωi 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 Di according to:
formula
(7)
where Ma and Pei 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:
formula
(8)
where the sign before the lattice-based Mach number is positive for forward (i.e. upward, k = 1) streaming or negative for backward (i.e. downward, k = 2) streaming.
With regard to inlet boundary conditions, one obtains via backward streaming from the adjacent lattice site while fi,1(0,t) is unknown. Danckwerts inlet condition is imposed by combining first-order finite-differences estimate of ∂yi/∂z together with Equation (4), thus leading to:
formula
(9)
where yi,in is the species concentration in the feeding effluent. Dirichlet inlet condition can be obtained by letting Pei →∞ in the above expression, namely:
formula
(10)
At exit, is obtained from forward streaming while fi,2(L,t) is unknown. Null Neumann condition is imposed by again combining first-order finite-differences estimate of ∂yi/∂z with Equation (4) so that the final result is:
formula
(11)
Last but not least, initial conditions for particle distribution functions are implemented as:
formula
(12)
where yi(z,0) are the initial conditions imposed to chemical species concentrations.
As initial development, the LBM simulator was implemented to numerically solve Equation (1) for species concentrations y1 to y5 (in Table 3) together with biomass concentration y7. In view of that, Table 4 shows corresponding generation/consumption rates as expressed in terms of a mutual auxiliary parameter A, namely:
formula
(13)
Table 4

Reaction rates, and initial and feeding concentrations for chemical species invoked in the LBM simulator

yiChemical speciesReaction rateInitial concentrationFeeding concentration
y1 COD  y1(z,0) = 0 y1,in = 24.7 kg-COD/m3 
y2 Acetic acid  y2(z,0) = 0 y2,in = 0.5013 kg-COD/m3 
y3 Propionic acid  y3(z,0) = 0 y3,in = 0.0801 kg-COD/m3 
y4 Butyric acid  y4(z,0) = 0 y4,in = 0.5266 kg-COD/m3 
y5 Dissolved H2  y5(z,0) = 0 y5,in = 0.0001 kg-COD/m3 
y7 Biomass  y7(z,0) = 0.1 kg-COD/m3 y7,in = 0.1 kg-COD/m3 
yiChemical speciesReaction rateInitial concentrationFeeding concentration
y1 COD  y1(z,0) = 0 y1,in = 24.7 kg-COD/m3 
y2 Acetic acid  y2(z,0) = 0 y2,in = 0.5013 kg-COD/m3 
y3 Propionic acid  y3(z,0) = 0 y3,in = 0.0801 kg-COD/m3 
y4 Butyric acid  y4(z,0) = 0 y4,in = 0.5266 kg-COD/m3 
y5 Dissolved H2  y5(z,0) = 0 y5,in = 0.0001 kg-COD/m3 
y7 Biomass  y7(z,0) = 0.1 kg-COD/m3 y7,in = 0.1 kg-COD/m3 

For the chemical species preliminary invoked, Table 4 also shows initial conditions yi(z,0) for Equation (12) as well as feeding values yi,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 kd1 = 0.02 day−1, Ysu = 0.1, fC2,su = 0.41, fC3,su = 0.27, fC4,su = 0.13, fH2,su = 0.19, and Di = 0.02 m2/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 y1 to y4 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.

Table 5

Comparing species concentrations measured at EC exit after 1 day of operation (Ferraz et al. 2014) with counterparts numerically simulated via LBM at z = L and t = 1 day

Species concentrationHRT = 24 h
HRT = 16 h
HRT = 12 h
HRT = 8 h
LBMMeasuredLBMMeasuredLBMMeasuredLBMMeasured
y1 (kg-COD/m38.8 19.6 ± 3.2 17.6 20.0 ± 1.7 17.3 20.9 ± 3.1 17.6 22.1 ± 2.2 
y2 (kg-COD/m32.1 2.0 ± 0.8 3.1 1.5 ± 0.7 3.2 1.7 ± 0.4 3.1 0.8 ± 0.2 
y3 (kg-COD/m31.2 0.8 ± 0.5 1.8 1.0 ± 0.3 1.9 1.1 ± 0.4 1.8 1.2 ± 0.1 
y4 (kg-COD/m30.9 0.9 ± 0.1 1.3 1.3 ± 0.5 1.4 1.5 ± 0.7 1.4 2.9 ± 0.6 
Species concentrationHRT = 24 h
HRT = 16 h
HRT = 12 h
HRT = 8 h
LBMMeasuredLBMMeasuredLBMMeasuredLBMMeasured
y1 (kg-COD/m38.8 19.6 ± 3.2 17.6 20.0 ± 1.7 17.3 20.9 ± 3.1 17.6 22.1 ± 2.2 
y2 (kg-COD/m32.1 2.0 ± 0.8 3.1 1.5 ± 0.7 3.2 1.7 ± 0.4 3.1 0.8 ± 0.2 
y3 (kg-COD/m31.2 0.8 ± 0.5 1.8 1.0 ± 0.3 1.9 1.1 ± 0.4 1.8 1.2 ± 0.1 
y4 (kg-COD/m30.9 0.9 ± 0.1 1.3 1.3 ± 0.5 1.4 1.5 ± 0.7 1.4 2.9 ± 0.6 

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 yi 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 yi(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 Di (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.

Figure 3

Time-dependent concentrations yi(L,t) simulated at APBR exit by imposing either Danckwerts (dots) or Dirichlet condition (solid line) at APBR inlet for chemical species invoked in the preliminary model: (a) COD (i = 1); (b) acetic acid (i = 2); (c) propionic acid (i = 3); (d) butyric acid (i = 4); (e) dissolved H2 (i = 5); (f) biomass (i = 7).

Figure 3

Time-dependent concentrations yi(L,t) simulated at APBR exit by imposing either Danckwerts (dots) or Dirichlet condition (solid line) at APBR inlet for chemical species invoked in the preliminary model: (a) COD (i = 1); (b) acetic acid (i = 2); (c) propionic acid (i = 3); (d) butyric acid (i = 4); (e) dissolved H2 (i = 5); (f) biomass (i = 7).

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 yi(z,t) numerically simulated over the solution domain (i.e. for 0 ≤ zL) at the following instants: t = 0.25 day (=tfinal/4), t = 0.5 day (=tfinal/2), t = 0.75 day (=3tfinal/4), and t = 1.0 day (=tfinal).

Figure 4

Species concentrations yi(z,t) throughout the solution domain at specific time instants simulated under Danckwerts condition at APBR inlet for chemical species invoked in the preliminary model: (a) COD (i = 1); (b) acetic acid (i = 2); (c) propionic acid (i = 3); (d) butyric acid (i = 4); (e) dissolved H2 (i = 5); (f) biomass (i = 7).

Figure 4

Species concentrations yi(z,t) throughout the solution domain at specific time instants simulated under Danckwerts condition at APBR inlet for chemical species invoked in the preliminary model: (a) COD (i = 1); (b) acetic acid (i = 2); (c) propionic acid (i = 3); (d) butyric acid (i = 4); (e) dissolved H2 (i = 5); (f) biomass (i = 7).

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 H2 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 y7(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 y7(0,t) remains fixed at the feeding value, i.e. y7,in = 0.1 kg-COD/m3 (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 Di) 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.

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/m3 for HRT = 24 h) and butyric acid (1.3 kg-COD/m3 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.

Bathe
,
K.-J.
1996
Finite Elements Procedures
.
Prentice-Hall
,
Upper Saddle River, NJ
,
USA
.
Batstone
,
D. J.
,
Keller
,
J.
,
Angelidaki
,
I.
,
Kalyuzhnyi
,
S. V.
,
Pavostathis
,
S. G.
,
Rozzi
,
A.
,
Sanders
,
W. T.
,
Siegrist
,
H.
&
Vavilin
,
V. A.
2002
The IWA Anaerobic Digestion Model No. 1 (ADM1)
.
Water Sci. Technol.
45
(
10
),
65
73
.
Brannock
,
M.
,
Leslie
,
G.
,
Wang
,
Y.
&
Buetehorn
,
S.
2010
Optimising mixing and nutrient removal in membrane bioreactors: CFD modelling and experimental validation
.
Desalination
250
(
2
),
815
818
.
Chen
,
X.
2010
Simulation of 2D cavitation bubble growth under shear flow by lattice Boltzmann model
.
Commun. Comput. Phys.
7
(
1
),
212
223
.
Datta
,
A. K.
,
Sablani
,
S. S.
2007
Mathematical modeling techniques in food and bioprocess: an overview
. In:
Handbook of Food and Bioprocess Modeling Techniques
(
Sablani
,
S. S.
,
Rahman
,
M. S.
,
Datta
,
A. K.
&
Mujumdar
,
A. R.
, eds).
CRC Press
,
Boca Raton, FL
,
USA
, pp.
1
11
.
de Souza-Santos
,
M. L.
2010
Solid Fuels Combustion and Gasification: Modeling, Simulation, and Equipment Operations
, 2nd edn.
CRC Press
,
Boca Raton, FL
,
USA
.
Durán
,
R.
,
Villa
,
A. L.
,
Ribeiro
,
R.
&
Rabi
,
J. A.
2015
Pectin extraction from mango peels in batch reactor: dynamic one-dimensional modeling and lattice Boltzmann simulation
.
Chem. Prod. Proc. Model.
10
(
3
),
203
210
.
Ferraz
,
A. D. N.
Jr
2013
Anaerobic Digestion of Sugar Cane Vinasse in Acidogenic Fixed Bed Reactor Followed by Methanogenic Reactor Sludge Blanket Type
.
PhD Thesis
(in Portuguese)
,
São Carlos School of Engineering, University of São Paulo
,
São Carlos
,
Brazil
.
Ferraz
,
A. D. N.
Jr
,
Wenzel
,
J.
,
Etchebehere
,
C.
&
Zaiat
,
M.
2014
Effect of organic loading rate on hydrogen production from sugarcane vinasse in thermophilic acidogenic packed bed reactors
.
Int. J. Hydrogen Energy
39
(
30
),
16852
16862
.
Gaspar
,
F.
,
Lu
,
T.
,
Santos
,
R.
&
Al-Duri
,
B.
2003
Modelling the extraction of essential oils with compressed carbon dioxide
.
J. Supercrit. Fluids
25
(
3
),
247
260
.
Głuszcz
,
P.
,
Petera
,
J.
&
Ledakowicz
,
S.
2011
Mathematical modeling of the integrated process of mercury bioremediation in the industrial bioreactor
.
Bioproc. Biosyst. Eng.
34
(
3
),
275
285
.
Kreyszig
,
E.
1993
Advanced Engineering Mathematics
. 7th edn.
John Wiley & Sons
,
New York
,
USA
.
McNamara
,
G. R.
&
Zanetti
,
G.
1988
Use of the Boltzmann equation to simulate lattice-gas automata
.
Phys. Rev. Lett.
61
(
20
),
2332
2335
.
Mohamad
,
A. A.
2011
Lattice Boltzmann Method: Fundamentals and Engineering Applications with Computer Codes
.
Springer-Verlag
,
London
,
UK
.
Naessens
,
W.
,
Maere
,
T.
,
Ratkovich
,
N. R.
,
Vedantam
,
S.
&
Nopens
,
I.
2012
Critical review of membrane bioreactor models: part 2: hydrodynamic and integrated models
.
Biores. Technol.
122
,
107
118
.
Parco
,
V.
,
Du Toit
,
G.
,
Wentzel
,
M.
&
Ekama
,
G.
2007
Biological nutrient removal in membrane bioreactors: denitrification and phosphorus removal kinetics
.
Water Sci. Technol.
56
(
6
),
125
134
.
Patankar
,
S. V.
1980
Numerical Heat Transfer and Fluid Flow
.
Hemisphere Publishing Corporation
,
Washington DC
,
USA
.
Picioreanu
,
C.
,
van Loosdrecht
,
M. C. M.
&
Heijnen
,
J. J.
2001
Two-dimensional model of biofilm detachment caused by internal stress from liquid flow
.
Biotechnol. Bioeng.
72
(
2
),
205
218
.
Qian
,
Y. H.
,
D'Humières
,
D.
&
Lallemand
,
P.
1992
Lattice BGK models for Navier-Stokes equation
.
Europhys. Lett.
17
,
479
484
.
Rabi
,
J. A.
2012
Numerical simulation of mass transfer aiming at food processes and bioprocesses in fixed-bed equipment: a comparison between finite-differences and lattice-Boltzmann method
.
WSEAS Trans. Heat Mass Transf.
7
,
69
78
.
Reis-Vasco
,
E. M. C.
,
Coelho
,
J. A. P.
,
Palavra
,
A. M. F.
,
Marrone
,
C.
&
Reverchon
,
E.
2000
Mathematical modelling and simulation of pennyroyal essential oil supercritical extraction
.
Chem. Eng. Sci.
55
(
15
),
2917
2922
.
Rosa
,
R. H.
,
von Atzingen
,
G. V.
,
Belandria
,
V.
,
Oliveira
,
A. L.
,
Bostyn
,
S.
&
Rabi
,
J. A.
2016
Lattice Boltzmann simulation of cafestol and kahweol extraction from green coffee beans in high-pressure system
.
J. Food Eng.
176
,
88
96
.
Sanchez
,
F. J. M.
,
del Valle
,
E. M.
,
Serrano
,
M. A. G.
&
Cerro
,
R. L.
2004
Modeling of monolith supported affinity chromatography
.
Biotechnol. Progr.
20
(
3
),
811
817
.
Štumbauer
,
V.
,
Petera
,
K.
&
Štys
,
D.
2013
The lattice Boltzmann method in bioreactor design and simulation
.
Math. Comp. Model.
57
(
7–8
),
1913
1918
.
Succi
,
S.
2001
The Lattice Boltzmann Equation for Fluid Dynamics and Beyond
.
Oxford University Press
,
New York
,
USA
.
Sukop
,
M. C.
&
Thorne
,
D. T.
Jr
2006
Lattice Boltzmann Modeling – An Introduction for Geoscientists and Engineers
.
Springer-Verlag
,
Berlin
,
Germany
.
van der Sman
,
R. G. M.
2007
Lattice Boltzmann simulation of microstructures
. In:
Handbook of Food and Bioprocess Modeling Techniques
(
Sablani
,
S. S.
,
Rahman
,
M. S.
,
Datta
,
A. K.
&
Mujumdar
,
A. R.
, eds).
CRC Press
,
Boca Raton, FL
,
USA
, pp.
15
39
.
Vlissidis
,
A.
&
Zouboulis
,
A. I.
1993
Thermophilic anaerobic digestion of alcohol distillery wastewaters
.
Bioresour. Technol.
43
(
2
),
131
140
.